Libraries

The following chunks are related to the article Green Bubbles: A Four-Stage Paradigm for Detection and Propagation. The project requires the libraries listed below (please ensure they are installed before running the code).

After installing all the required packages, the libraries are loaded.

library(bbdetection)
library(bruceR)
## 
## bruceR (v2025.8)
## Broadly Useful Convenient and Efficient R functions
## 
## Packages also loaded:
## ✔ dplyr      ✔ data.table
## ✔ tidyr      ✔ emmeans
## ✔ stringr    ✔ lmerTest
## ✔ forcats    ✔ effectsize
## ✔ ggplot2    ✔ performance
## ✔ cowplot    ✔ interactions
## 
## Main functions of `bruceR`:
## cc()             Describe()  TTEST()
## add()            Freq()      MANOVA()
## .mean()          Corr()      EMMEANS()
## set.wd()         Alpha()     PROCESS()
## import()         EFA()       model_summary()
## print_table()    CFA()       lavaan_summary()
## 
## For full functionality, please install all dependencies:
## install.packages("bruceR", dep=TRUE)
## 
## Online documentation:
## https://psychbruce.github.io/bruceR
## 
## To use this package in publications, please cite:
## Bao, H. W. S. (2021). bruceR: Broadly useful convenient and efficient R functions (Version 2025.8) [Computer software]. https://doi.org/10.32614/CRAN.package.bruceR
## 
## These packages are dependencies but not yet installed:
## - pacman, openxlsx, ggtext, see, phia, MuMIn, BayesFactor, GGally
## 
## ***** Install All Dependencies *****
## install.packages("bruceR", dep=TRUE)
library(bsts)
## Loading required package: BoomSpikeSlab
## Loading required package: Boom
## 
## Attaching package: 'Boom'
## The following object is masked from 'package:stats':
## 
##     rWishart
## 
## Attaching package: 'BoomSpikeSlab'
## The following object is masked from 'package:stats':
## 
##     knots
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:data.table':
## 
##     yearmon, yearqtr
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: xts
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:data.table':
## 
##     first, last
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## Attaching package: 'bsts'
## The following object is masked from 'package:BoomSpikeSlab':
## 
##     SuggestBurn
library(cpm)
library(dplyr)
library(exuber)
## 
## Attaching package: 'exuber'
## The following objects are masked from 'package:zoo':
## 
##     index, index<-
library(fDMA)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## 
## Attaching package: 'fDMA'
## The following object is masked from 'package:effectsize':
## 
##     standardize
library(forecast)
library(ggfortify)
## Registered S3 methods overwritten by 'ggfortify':
##   method                 from    
##   autoplot.Arima         forecast
##   autoplot.acf           forecast
##   autoplot.ar            forecast
##   autoplot.bats          forecast
##   autoplot.decomposed.ts forecast
##   autoplot.ets           forecast
##   autoplot.forecast      forecast
##   autoplot.stl           forecast
##   autoplot.ts            forecast
##   fitted.ar              forecast
##   fortify.ts             forecast
##   residuals.ar           forecast
library(ggplot2)
library(ghyp)
## Loading required package: numDeriv
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(grid)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:cowplot':
## 
##     stamp
## The following objects are masked from 'package:data.table':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(MCS)
library(Metrics)
## 
## Attaching package: 'Metrics'
## The following object is masked from 'package:forecast':
## 
##     accuracy
## The following objects are masked from 'package:performance':
## 
##     mae, mse, rmse
library(multDM)
library(readxl)
library(rugarch)
## Loading required package: parallel
## 
## Attaching package: 'rugarch'
## The following object is masked from 'package:exuber':
## 
##     report
library(tseries)
library(urca)
library(vars)
## Loading required package: strucchange
## Loading required package: sandwich
## 
## Attaching package: 'strucchange'
## The following object is masked from 'package:stringr':
## 
##     boundary
## Loading required package: lmtest
library(xts)
library(zoo)

Data generation Process and simulations

Firstly, the data generation of bubbles phenomena is outlined (this step may take a few seconds).

Three different bubble patterns are being considered.

Below is a graphical representation of potential patterns of bubble collapse as reported in Phillips and Shi (2018).

The Backward Sup ADF(BSADF) is compared with the Kolmogorov-Smirnov Change Point Detection Model (KS-CPM). The disturbing path is initially examined, followed by 1000 simulations (this process may take a few seconds).

Below, the simulation results are summarized, taking into account both the frequency of correct bubble identifications and the Root Mean Square Error (RMSE). This is presented in Table A.1 (a) of the article.

#Correctness KS-CPM
matrix_test = matrix(ncol = 3, nrow = correct)
matrix_ADF = matrix(ncol = 3, nrow = correct_ADF)
for (i in 1:correct) {
  matrix_test[i, ] <-
    c(results_DIS[lengths(results_DIS) == 3][[i]][1], results_DIS[lengths(results_DIS) == 3] [[i]][2], results_DIS[lengths(results_DIS) == 3][[i]][3])
}
cat("Correctness KS-CPM:", dim(matrix_test)[1], "\n")
## Correctness KS-CPM: 906
#Correctness BSADF test
for (i in 1:correct_ADF) {
  matrix_ADF[i, ] <-
    c(as.numeric(substr(ADF_DIS[lengths(ADF_DIS) == 3][[i]][1], 1, 2)),
      as.numeric(substr(ADF_DIS[lengths(ADF_DIS) == 3][[i]][2], 1, 2)),
      as.numeric(substr(ADF_DIS[lengths(ADF_DIS) == 3][[i]][3], 1, 2)))
}
cat("Correctness BSADF:", dim(matrix_ADF)[1], "\n")
## Correctness BSADF: 854
#RMSE KS-CPM
cat("RMSE KS-CPM:", round(Metrics::rmse(
  c(matrix_test[, 1], matrix_test[, 2], matrix_test[, 3]),
  c(rep(40, correct), rep(60, correct), rep(70, correct))
), 3), "\n")
## RMSE KS-CPM: 1.884
#RMSE BSADF test
cat("RMSE BSADF:", round(Metrics::rmse(
  c(matrix_ADF[, 1], matrix_ADF[, 2], matrix_ADF[, 3]), c(
    rep(40, correct_ADF),
    rep(60, correct_ADF),
    rep(70, correct_ADF)
  )
), 3), "\n")
## RMSE BSADF: 4.531

The same analysis is conducted for the abrupt or sudden path (this may take a few seconds, with 1,000 replications carried out).

This is presented in Table A.1(b) of the article.

#Correctness KS-CPM
matrix_test = matrix(ncol = 3, nrow = correct)
matrix_ADF = matrix(ncol = 3, nrow = correct_ADF)
for (i in 1:correct) {
  matrix_test[i, ] <-
    c(results_SUD[lengths(results_SUD) == 3][[i]][1], results_SUD[lengths(results_SUD) == 3] [[i]][2], results_SUD[lengths(results_SUD) == 3] [[i]][3])
}
cat("Correctness KS-CPM:", dim(matrix_test)[1], "\n")
## Correctness KS-CPM: 879
#Correctness BSADF test
for (i in 1:correct_ADF) {
  matrix_ADF[i, ] <-
    c(as.numeric(substr(ADF_SUD[lengths(ADF_SUD) == 3][[i]][1], 1, 2)),
      as.numeric(substr(ADF_SUD[lengths(ADF_SUD) == 3][[i]][2], 1, 2)),
      as.numeric(substr(ADF_SUD[lengths(ADF_SUD) == 3][[i]][3], 1, 2)))
}
cat("Correctness BSADF:", dim(matrix_ADF)[1], "\n")
## Correctness BSADF: 566
#RMSE KS-CPM
cat("RMSE KS-CPM:", round(Metrics::rmse(
  c(matrix_test[, 1], matrix_test[, 2], matrix_test[, 3]),
  c(rep(50, correct), rep(60, correct), rep(66, correct))
), 3), "\n")
## RMSE KS-CPM: 1.912
#RMSE BSADF test
cat("RMSE BSADF:", round(Metrics::rmse(
  c(matrix_ADF[, 1], matrix_ADF[, 2], matrix_ADF[, 3]), c(
    rep(50, correct_ADF),
    rep(60, correct_ADF),
    rep(66, correct_ADF)
  )
), 3), "\n")
## RMSE BSADF: 3.119

Finally, the analysis is extended to include the smoothing path as well (this may take a few seconds, with 1,000 replications carried out).

This is presented in Table A.1(c) of the article.

#Correctness KS-CPM
matrix_test = matrix(ncol = 3, nrow = correct)
matrix_ADF = matrix(ncol = 3, nrow = correct_ADF)
for (i in 1:correct) {
  matrix_test[i, ] <-
    c(results_SMO[lengths(results_SMO) == 3][[i]][1], results_SMO[lengths(results_SMO) == 3] [[i]][2], results_SMO[lengths(results_SMO) == 3] [[i]][3])
}
cat("Correctness KS-CPM:", dim(matrix_test)[1], "\n")
## Correctness KS-CPM: 892
#Correctness BSADF test
for (i in 1:correct_ADF) {
  matrix_ADF[i, ] <-
    c(as.numeric(substr(ADF_SMO[lengths(ADF_SMO) == 3][[i]][1], 1, 2)),
      as.numeric(substr(ADF_SMO[lengths(ADF_SMO) == 3][[i]][2], 1, 2)),
      as.numeric(substr(ADF_SMO[lengths(ADF_SMO) == 3][[i]][3], 1, 2)))
}
cat("Correctness BSADF:", dim(matrix_ADF)[1], "\n")
## Correctness BSADF: 812
#RMSE KS-CPM
cat("RMSE KS-CPM:", round(Metrics::rmse(
  c(matrix_test[, 1], matrix_test[, 2], matrix_test[, 3]),
  c(rep(40, correct), rep(60, correct), rep(80, correct))
), 3), "\n")
## RMSE KS-CPM: 2.223
#RMSE BSADF test
cat("RMSE BSADF:", round(Metrics::rmse(
  c(matrix_ADF[, 1], matrix_ADF[, 2], matrix_ADF[, 3]), c(
    rep(40, correct_ADF),
    rep(60, correct_ADF),
    rep(80, correct_ADF)
  )
), 3), "\n")
## RMSE BSADF: 6.138

The simulation study provides evidence that the change point detection model outperforms the BSADF test. With these results in hand, attention can now turn to the analysis of real-world data.

Detection phase

For the empirical analysis, the Renixx index is used as a real-world data example.

#NOTE: These plots are for reference only and not included in the article.
#Updating of the data.
renixx <- read.csv("Economic variables/Renixx.csv")
#Descriptive statistics
renixx$date <- ymd(renixx$date)
renixx <- renixx[order(renixx$date), ]
rownames(renixx) <- NULL
renixx <- renixx[renixx$date < '2024-03-01' &
                   renixx$date > '2004-12-31', ]
ggplot(renixx, aes(date, close)) + geom_line() + ggtitle('') +
  xlab('Time') + ylab('Closed Value')

#Log-returns
df <- renixx %>% dplyr::select(1, 3)
rownames(df) <- NULL
df$abs <- abs(c(0, diff(log(df$close))))
#Absolute log-returns
df$rate <- c(0, diff(log(df$close)))
ggplot(df, aes(date, abs)) + geom_line() + geom_line() + ggtitle('') +
  xlab('Time') + ylab('Absolute Returns')

The empirical analysis is conducted at a monthly frequency, as the forecasting phase incorporates textual data (e.g. Google Trends), which are available only at that frequency.

#NOTE: These plots are for reference only and not included in the article.
Month <- as.Date(cut(df$date, "month"))
renixx_m <- aggregate(close ~ Month, df, mean)
#Log-returns
renixx_m$Renixx <- c(0, diff(log(renixx_m$close)))
#Absolute log-returns
renixx_m$abs <- c(0, abs(diff(log(renixx_m$close))))
#Standardized values
renixx_m$St <- log(renixx_m$close)
ggplot(renixx_m, aes(Month, abs)) + geom_line() + xlab('Time') + ylab('Absolute Returns')

ggplot(renixx_m, aes(Month, close)) + geom_line() + xlab('Time') + ylab('Closed Value')

#NOTE: This plot is for reference only and not included in the article.
#KS-CPM
test <-
  processStream(
    ts(renixx_m$Renixx, freq = 12, start = decimal_date(ymd("2005-1-1"))),
    "Kolmogorov-Smirnov",
    ARL0 = 500,
    startup = 46
  ) #Startup = 46 (20%) and ARL0 = 500 are default values.
breack <- renixx_m$Month[test$changePoints]
plot(
  renixx_m$Month,
  renixx_m$close,
  type = "l",
  xlab = "Time",
  ylab = "Closed Value",
  bty = "l",
  main = ''
)
abline(v = breack, lty = 2, col = 'red')

test
## $x
##                Jan           Feb           Mar           Apr           May
## 2005  0.0000000000  0.0530189635  0.0687943955 -0.0119283046 -0.0031614423
## 2006  0.1576983710  0.1749323966  0.0341515622  0.0957208036 -0.0217993066
## 2007  0.0545578612  0.1358327830  0.0024174940  0.0961005547  0.0151803929
## 2008 -0.1756064749 -0.0968326005 -0.0728707695  0.1489470465  0.0631577298
## 2009  0.0914962218 -0.0611492288 -0.1513912861  0.1802162711  0.1857942118
## 2010 -0.0009596813 -0.1178961944 -0.0094209776  0.0579208666 -0.1142693067
## 2011  0.0444370211  0.0199152329 -0.0027564990 -0.0046400509 -0.0842015051
## 2012  0.0668437994  0.0400594438 -0.1211641718 -0.1309569248 -0.1264133545
## 2013  0.1169760667  0.0208876124  0.0235850195  0.0289741272  0.1650370226
## 2014  0.1084627881  0.0490217620  0.0496114001 -0.0393263131  0.0018409527
## 2015  0.0505690299  0.1119845016  0.1203081711  0.0974680083  0.0076100111
## 2016 -0.0847890610 -0.0993756016  0.0805477489  0.0068269693 -0.0454210495
## 2017  0.0467343606  0.0335219319 -0.0036486303  0.0306472072  0.0116983857
## 2018  0.0308244477 -0.0387336697  0.0138115033  0.0121550320  0.0402637115
## 2019  0.0212828655  0.0875994637  0.0182882461  0.0331691402 -0.0172071968
## 2020  0.0950428800  0.1414182560 -0.2384290135  0.0397472112  0.1202455451
## 2021  0.2610942237 -0.0449369726 -0.1969225402 -0.0142498162 -0.1437374275
## 2022 -0.1474325870 -0.0577322742  0.1585923453  0.0035731285 -0.1323672376
## 2023 -0.0219598274  0.0442065571 -0.0351862986 -0.0386876920 -0.0410413230
## 2024 -0.0509250599 -0.0444162248                                          
##                Jun           Jul           Aug           Sep           Oct
## 2005  0.1247902752  0.0454159579  0.0543622625  0.1682175365 -0.0217539966
## 2006 -0.1652615193 -0.0083861601  0.0142407772 -0.0494610158  0.0040148755
## 2007  0.0444541078  0.0909225030 -0.0643837096  0.0545671136  0.1226449466
## 2008 -0.0137677702 -0.0973702759  0.0520610621 -0.0854028706 -0.4664049460
## 2009 -0.0332124142 -0.0552182147 -0.0052270289  0.0006166750 -0.0346679902
## 2010 -0.0227131829  0.0187390404 -0.0285953805 -0.0115045088 -0.0298951698
## 2011 -0.0882440279 -0.0532377107 -0.2019251835 -0.1352156857 -0.1116346573
## 2012 -0.1054888691 -0.0305734416  0.0399615938  0.0434881288 -0.0568106746
## 2013 -0.0042486901  0.0675520572  0.0277323545  0.0561515696  0.1365425298
## 2014  0.0789394331  0.0028426798  0.0254331442  0.0262621995 -0.1329449568
## 2015 -0.0156086339 -0.0041899030 -0.1000522770 -0.0794318098  0.0401108001
## 2016 -0.0051328820  0.0255196141  0.0241075768 -0.0108726803  0.0041885525
## 2017  0.0126864317 -0.0176751480 -0.0165971458  0.0071993065  0.0379319717
## 2018 -0.0562332544 -0.0113617454  0.0067924042 -0.0288366898 -0.0402676312
## 2019  0.0260556631  0.0598687114  0.0045555167  0.0363606959 -0.0027454836
## 2020  0.0796221242  0.1632228767  0.1336612166  0.0380903691  0.1797831938
## 2021  0.0706035278  0.0517403319  0.0110938875  0.0005457772  0.0582821748
## 2022  0.0420121928  0.0639393711  0.1510654697 -0.0462191791 -0.1653198471
## 2023  0.0391340737 -0.0103654819 -0.1302217746 -0.0716099572 -0.1255400396
## 2024                                                                      
##                Nov           Dec
## 2005  0.0008430376 -0.0362914995
## 2006  0.0806226741  0.0677080589
## 2007  0.0690391310  0.0735464512
## 2008 -0.2320893945 -0.0344547376
## 2009 -0.0540133749  0.0410858613
## 2010 -0.0711144351 -0.0051122271
## 2011 -0.0953715468 -0.0654085981
## 2012 -0.0718484919  0.0734265032
## 2013  0.0610456880 -0.0599608792
## 2014  0.0379275901 -0.0498675522
## 2015 -0.0060524177  0.0324712978
## 2016 -0.0762955978  0.0166432637
## 2017 -0.0169453040 -0.0037327366
## 2018  0.0556210590 -0.0101469641
## 2019  0.0452281751  0.0595370791
## 2020  0.0707716327  0.1290430288
## 2021  0.1227725546 -0.1391067224
## 2022  0.0582839843 -0.0243450236
## 2023 -0.0499224818  0.1019640021
## 2024                            
## 
## $changePoints
## [1]  36  95 124 178 193
## 
## $detectionTimes
## [1]  71 101 161 182 226

A robustness analysis is subsequently conducted using weekly frequency data to confirm the results obtained from the low-frequency analysis.

#NOTE: This plot is for reference only and not included in the article
#KS-CPM
options(xts.message.period.apply.mean = FALSE)
weekly_dates <- floor_date(renixx$date, unit = "week", week_start = 1)
renixx_w <- apply.weekly(renixx, mean)
renixx_w$date <- unique(weekly_dates)
test_w <-
  processStream(c(0, diff(log(renixx_w$close))),
                "Kolmogorov-Smirnov",
                ARL0 = 5000,
                startup = 200) #Startup = 200 (20%) and ARL0 = 5000 are default values.
breack_w <- renixx_w$date[test_w$changePoints]
plot(
  renixx_w$date,
  renixx_w$close,
  type = "l",
  xlab = "Time",
  ylab = "Closed Value",
  bty = "l",
  main = ''
)
abline(v = breack_w, lty = 2, col = 'red')

test_w
## $x
##    [1]  0.000000e+00 -5.878894e-03  1.953790e-02 -4.480740e-03  1.966467e-02
##    [6]  1.918462e-02  1.270518e-02  1.400954e-02  5.259928e-02  2.431125e-03
##   [11] -7.395833e-03  5.205206e-03 -1.669158e-03  3.259430e-02 -2.238574e-02
##   [16] -5.804813e-02 -2.877861e-03  3.948530e-03  1.598808e-02  2.247469e-02
##   [21]  1.878275e-02  3.566613e-02  3.030427e-02  4.466098e-02  2.633549e-02
##   [26] -1.677556e-02  1.721606e-02  1.522883e-02 -1.932343e-02  2.670012e-02
##   [31]  2.876340e-02 -8.398668e-03  1.249381e-02  1.071722e-02  7.129308e-02
##   [36]  5.944531e-02  2.904082e-02  2.375016e-02  2.094917e-02  1.089909e-02
##   [41] -5.708736e-02 -5.482698e-02 -8.122804e-03  2.962571e-02  2.120511e-02
##   [46]  1.424928e-02  3.541012e-04 -6.393643e-02 -6.633199e-03  7.751800e-03
##   [51]  1.939477e-02  1.369966e-02  3.426673e-02  2.793011e-02  9.378945e-02
##   [56]  5.177254e-02  7.836332e-02  4.082264e-03 -4.174377e-03  4.316745e-02
##   [61]  3.452686e-02 -2.984601e-02  3.113674e-03 -1.131846e-02  3.304017e-02
##   [66]  3.405104e-02  1.184320e-02  7.252018e-02  3.189595e-04  2.273419e-02
##   [71]  2.704319e-02 -1.458327e-01 -7.410913e-02  2.016918e-02 -5.634466e-02
##   [76] -7.035113e-02  4.419808e-02  1.411004e-02  5.626760e-02 -6.753917e-02
##   [81] -4.827715e-02  1.115619e-02  3.707433e-02  4.275845e-03  1.209622e-02
##   [86] -1.255529e-04 -1.163105e-02 -1.279424e-02 -2.811397e-02 -2.482298e-02
##   [91] -6.719437e-03 -6.519007e-03  1.849707e-02  3.779415e-02  1.462364e-02
##   [96] -2.385967e-02  3.596141e-02  4.100056e-02  3.342278e-02  8.368814e-03
##  [101]  6.440725e-03  1.063366e-02  1.514957e-02  2.756914e-03  1.898669e-03
##  [106]  4.246108e-03  2.953157e-02  4.192369e-02  2.781282e-02  5.882827e-02
##  [111]  2.198549e-02  1.993698e-02 -5.644978e-02 -3.420726e-03  1.947554e-02
##  [116]  3.607947e-02  1.615364e-02  1.740880e-02  3.077813e-02  2.549127e-02
##  [121]  9.462664e-04 -1.710113e-02  4.369701e-03  2.232149e-03  1.958839e-02
##  [126]  1.929402e-02 -7.045796e-03  1.005621e-02  2.721616e-02 -1.398308e-03
##  [131]  6.327505e-02  3.056950e-02 -4.619331e-03 -2.670879e-02 -1.572936e-02
##  [136]  4.613523e-03 -7.633823e-02  1.697521e-02  1.558970e-02  4.658265e-02
##  [141] -3.038501e-02  3.302088e-02  6.856182e-02  2.473380e-02  2.291188e-02
##  [146]  6.327408e-03  2.146747e-02  6.809952e-02  5.735111e-02 -5.519734e-02
##  [151] -9.264226e-02  6.422694e-02  5.323074e-02  4.734835e-02 -9.395933e-03
##  [156]  4.415566e-02 -1.520653e-02 -1.054283e-01 -1.288330e-01 -1.290267e-01
##  [161]  5.998886e-02 -1.195138e-02  2.946524e-02 -2.236532e-03 -4.277900e-02
##  [166] -3.895420e-02 -1.662050e-02 -4.031452e-02  8.378635e-02  6.486860e-02
##  [171]  3.468499e-02  1.684634e-02  2.431819e-02  2.802520e-03 -3.407609e-03
##  [176]  5.162620e-02  1.854794e-02 -3.422227e-02  6.241718e-03 -4.761153e-02
##  [181]  5.413623e-02 -2.879984e-02 -9.029724e-02 -1.010833e-02  4.395118e-03
##  [186]  7.147566e-03  1.131879e-02 -8.585298e-03  4.867557e-03  6.486621e-02
##  [191]  3.715547e-02 -4.026662e-02 -1.085220e-01 -4.559488e-02  3.840129e-02
##  [196] -1.424247e-01 -2.725077e-01 -9.563063e-03 -1.131284e-01 -1.754655e-01
##  [201]  1.897018e-01 -1.608786e-01 -1.536952e-01  9.758902e-02 -2.561307e-03
##  [206]  1.273235e-02  3.274017e-03  2.656671e-02  4.767386e-02  9.252748e-02
##  [211] -8.880718e-02 -2.770057e-02  1.932301e-02 -6.300535e-03  2.271493e-02
##  [216] -4.507463e-02 -1.060467e-01 -9.234539e-02  2.685966e-03  1.834671e-02
##  [221]  7.544289e-02  4.564262e-02  2.824779e-02  6.832760e-02 -1.514913e-02
##  [226]  6.034581e-02  1.433860e-01 -4.762541e-02  2.866448e-02 -4.114632e-03
##  [231]  4.404866e-04  4.112824e-03 -4.077323e-02 -4.910487e-02  1.482071e-02
##  [236] -6.203254e-02  3.305877e-03  6.068155e-02  3.375006e-02 -1.991544e-02
##  [241] -9.662346e-03 -5.238441e-02 -8.562675e-03 -3.757045e-02  6.305983e-02
##  [246]  3.169411e-02  5.884604e-03 -2.885119e-02 -1.503463e-02 -5.851589e-03
##  [251] -1.700956e-02 -3.896213e-02 -3.323175e-02  2.396686e-02  1.586491e-02
##  [256] -2.051985e-02  9.627047e-03  3.341759e-02  7.918949e-03 -3.604328e-03
##  [261]  8.614766e-03  4.269299e-02 -1.928364e-02 -5.955869e-02 -5.913448e-02
##  [266] -7.185999e-03 -4.056594e-02  2.004755e-02 -4.066575e-02  2.040538e-03
##  [271]  1.813668e-02 -2.990231e-03  1.229888e-02  2.008199e-02  2.831490e-02
##  [276]  1.356334e-02 -1.193349e-02 -5.748650e-03 -5.215041e-02 -2.299303e-02
##  [281] -6.536692e-02 -2.417847e-02  1.765011e-02 -1.976185e-02  6.719968e-02
##  [286] -1.422085e-03 -6.497999e-02  3.249329e-02  2.622818e-02  1.234633e-02
##  [291]  5.910368e-03 -1.850083e-02 -6.971287e-03 -2.756235e-02 -4.426032e-02
##  [296]  2.139978e-02  2.244933e-02  8.491768e-04 -1.949946e-02 -8.471265e-03
##  [301] -9.756650e-03  9.264661e-03 -1.347088e-02 -1.645759e-02 -2.526974e-02
##  [306]  1.348098e-02 -5.346681e-02 -2.651538e-02 -7.929938e-03  3.333551e-02
##  [311]  2.038509e-02  1.200393e-03 -1.120975e-02  2.282231e-02  2.396949e-02
##  [316]  3.758950e-03 -2.314849e-03  3.104410e-03  5.451507e-04  3.018541e-02
##  [321] -2.370943e-02 -2.861251e-02 -2.533360e-02  6.350803e-02  3.247991e-03
##  [326]  5.603428e-02 -2.502147e-02 -3.951440e-02 -4.303812e-02 -6.065476e-03
##  [331] -2.430830e-02  4.923072e-04 -1.496862e-02 -4.021671e-02 -1.133597e-03
##  [336] -4.187488e-02 -9.983036e-03 -2.491424e-02  1.296504e-02  1.135810e-02
##  [341] -5.052253e-02 -2.070943e-02 -3.109689e-02 -6.605446e-02 -1.100289e-01
##  [346]  9.696869e-03 -3.746211e-02  4.123374e-02 -6.336236e-02 -3.156313e-02
##  [351] -7.999648e-02 -4.712305e-02 -3.367170e-02  2.062776e-02 -1.364994e-02
##  [356]  4.133833e-02 -4.918285e-02 -1.567954e-02 -4.819563e-02 -1.058920e-01
##  [361]  4.412089e-02  2.733836e-02 -5.577950e-02 -1.996030e-02  1.956732e-02
##  [366]  2.867973e-02  5.193022e-02  2.602178e-02 -1.029874e-02  2.350895e-02
##  [371]  4.084633e-02 -3.971194e-02 -9.754649e-03 -6.277963e-02 -5.075290e-02
##  [376]  2.427177e-02 -2.677284e-02 -2.741438e-02 -5.304858e-02 -3.294808e-02
##  [381] -9.716080e-03 -4.721272e-02 -6.016348e-03 -2.516360e-02 -5.502196e-02
##  [386] -5.932080e-02 -1.036819e-02 -3.780617e-02 -6.661474e-03  3.142222e-02
##  [391] -2.074170e-02  2.866192e-02 -4.103612e-02 -1.746855e-02 -5.856902e-02
##  [396]  1.799571e-02  6.881278e-02  1.070389e-02  3.374001e-02 -4.487528e-03
##  [401]  1.854513e-04  1.813460e-02  1.022491e-02 -1.130672e-02 -3.067019e-02
##  [406] -1.349574e-02 -1.728176e-03 -3.874606e-02 -1.866148e-02  8.278547e-03
##  [411] -3.779237e-02 -4.216367e-02  4.569670e-02  1.700872e-02  4.760759e-02
##  [416]  3.648354e-02  2.142344e-03  3.011089e-02  7.719132e-02 -3.029859e-03
##  [421]  1.253034e-03 -4.469717e-02  7.858996e-03  5.357637e-02  9.823990e-03
##  [426] -1.652110e-02  8.062175e-03  1.148753e-02  4.729153e-03 -1.100114e-02
##  [431] -2.371875e-02  4.685779e-02 -3.310412e-03  5.177142e-02  3.352090e-02
##  [436]  2.700745e-02  7.239754e-02  3.295036e-02  2.197115e-02 -3.864856e-02
##  [441] -2.388485e-02 -1.415738e-02 -1.576902e-02  5.094656e-02  2.985005e-02
##  [446]  1.945296e-02  1.546956e-02  1.518047e-02  4.394072e-03 -2.868809e-03
##  [451] -2.237068e-02 -1.294166e-02  3.961627e-02  3.414711e-02  1.282358e-02
##  [456]  1.756735e-02  5.757328e-02  2.056501e-02  5.046731e-02  2.923628e-02
##  [461] -3.182924e-02  5.763945e-02  1.668067e-03 -1.673510e-02 -1.013587e-02
##  [466] -2.447241e-02 -4.285748e-02  1.298133e-02  2.243867e-02  3.457966e-02
##  [471]  4.881900e-02  2.935455e-02  4.960184e-03 -3.541674e-02  2.294954e-03
##  [476]  5.175799e-02  1.970118e-02  3.905385e-02  2.929175e-02 -2.607381e-02
##  [481] -5.174362e-03 -3.185605e-02  2.836734e-02 -2.684288e-02 -4.176563e-02
##  [486]  3.945670e-02 -2.791146e-02 -1.593267e-03  1.292533e-02  1.111383e-02
##  [491]  4.392791e-02  4.771248e-03  1.798719e-02  2.956714e-02  3.163882e-03
##  [496]  1.458436e-02 -4.049437e-02 -5.912530e-03  1.310336e-02 -3.954271e-03
##  [501] -1.325636e-02  3.994567e-02  2.575020e-02  1.705765e-02  1.433792e-02
##  [506]  2.228403e-03 -4.069657e-02 -1.634728e-02 -2.600290e-02 -5.492875e-02
##  [511] -8.681718e-02  6.178665e-02  1.549733e-02  3.095362e-02 -2.100161e-02
##  [516]  5.517047e-03  2.411391e-02 -4.675571e-02 -2.005972e-02 -2.465975e-02
##  [521]  2.943684e-02  1.913005e-02  7.366457e-03  2.755054e-03  7.151800e-03
##  [526]  6.210496e-02  3.173737e-02  2.061677e-02  1.176297e-02  2.944299e-02
##  [531]  5.487867e-02  3.031556e-02  1.994988e-02 -1.250837e-02  3.049927e-02
##  [536]  5.354909e-02  1.957533e-02  1.095694e-03 -8.006187e-03 -2.655408e-02
##  [541]  2.476649e-02  1.518145e-02  1.928028e-02 -3.022108e-02 -2.726844e-02
##  [546]  2.826975e-03  2.470292e-02 -2.891986e-02 -2.860660e-02  5.983606e-02
##  [551]  8.151155e-03 -4.184714e-02 -4.878848e-03 -4.855978e-02 -4.320754e-02
##  [556] -9.787628e-02  2.818549e-02  1.718006e-02 -2.738252e-03 -1.950689e-02
##  [561] -3.247750e-02  8.456847e-02 -4.105144e-03 -3.155711e-03  9.689400e-03
##  [566]  6.974838e-03 -2.790823e-02 -1.616001e-02  1.868171e-02  3.091387e-02
##  [571] -3.799421e-02  3.945657e-02  3.705684e-02 -3.923505e-04 -3.068984e-02
##  [576] -9.367308e-02 -5.436506e-02  2.845817e-02 -2.089512e-02 -1.099003e-01
##  [581]  4.352881e-02  3.842879e-02  4.336906e-02  9.987669e-03  1.789261e-02
##  [586] -7.557119e-04 -8.816401e-03 -8.170221e-03  1.336398e-02  9.536530e-03
##  [591] -1.607636e-03 -5.333351e-02 -1.629295e-02  5.069556e-03  2.757479e-02
##  [596]  1.527593e-02 -5.752544e-03 -5.447050e-02  8.545192e-03  2.948988e-03
##  [601]  1.163174e-02  2.473625e-02  1.273662e-02  9.766889e-03  1.236243e-03
##  [606]  2.713442e-03 -5.043128e-03  1.168657e-02 -1.419903e-02  9.046009e-03
##  [611] -2.800544e-02  1.367153e-02  7.232231e-03  2.549408e-03 -1.752511e-02
##  [616]  2.007730e-02 -1.154641e-03 -5.278231e-02 -3.610572e-02 -1.334449e-02
##  [621]  1.998077e-02  9.625630e-03 -1.207688e-02  1.263553e-02  1.903079e-02
##  [626]  1.591883e-02  1.613147e-02  2.999510e-03  3.546153e-03 -2.984155e-03
##  [631] -2.216767e-03  1.572729e-02  3.374876e-02  6.893422e-03 -1.619556e-02
##  [636] -1.307154e-02  1.504801e-03 -6.088355e-03  1.660909e-02  2.762179e-02
##  [641]  7.131814e-03 -2.063085e-02  5.393492e-04  1.108007e-02  2.269518e-02
##  [646] -1.866670e-02  6.137448e-04  5.976420e-03  9.041692e-03  1.472981e-03
##  [651]  8.155439e-03 -1.006334e-02 -2.011606e-02  1.476951e-03  1.048418e-02
##  [656] -4.100339e-03 -2.184191e-02  1.598091e-02 -7.286760e-03 -1.519224e-02
##  [661] -3.353584e-03  3.407947e-03  2.011853e-02  4.067568e-03 -4.298471e-03
##  [666]  2.266538e-02  1.065335e-02 -5.886460e-04  2.198133e-03  1.743855e-02
##  [671] -1.419740e-02 -2.550876e-02 -1.502153e-02 -5.011513e-03 -8.292950e-03
##  [676]  3.534811e-02  1.012910e-02 -1.226374e-03  3.896401e-03  2.163656e-02
##  [681]  7.906036e-04 -1.725274e-02 -1.104366e-02 -4.657273e-02  1.111470e-02
##  [686]  2.823634e-02  6.953545e-03 -2.139619e-02  2.701648e-02 -2.925272e-04
##  [691] -2.521268e-02  1.089412e-02  1.260901e-02  2.500900e-03  3.039138e-03
##  [696]  2.102341e-02  1.661929e-02  1.084638e-02 -5.802673e-03 -1.990186e-02
##  [701] -3.674694e-02 -8.783411e-04 -4.053830e-03 -2.390901e-02 -3.540318e-03
##  [706]  1.800095e-02  1.947000e-03  6.160779e-05  1.635618e-03  1.174850e-02
##  [711] -1.801527e-02  2.185928e-03 -2.345092e-03 -2.049015e-02 -8.834913e-03
##  [716]  7.023527e-03  3.182423e-03  1.342155e-03 -5.790531e-02  1.648774e-03
##  [721] -1.658312e-03  7.568649e-03  5.692725e-02  1.393352e-02 -2.043806e-02
##  [726]  3.037158e-02  9.835817e-03 -1.402710e-02 -3.707668e-02 -3.778128e-02
##  [731]  7.250161e-03  4.691653e-02  1.860220e-02  1.608158e-02  3.741827e-02
##  [736]  1.967001e-02  8.311829e-03  8.705086e-03  9.521085e-03 -3.325868e-03
##  [741]  1.581319e-02 -2.922110e-03 -1.352165e-02  3.548069e-02  4.628357e-03
##  [746]  1.349104e-03  4.261465e-03 -4.065192e-03 -8.063517e-03 -1.205098e-02
##  [751] -4.994066e-03 -1.725314e-03  7.551257e-03  2.749150e-02  1.404949e-02
##  [756] -2.711345e-03  3.169706e-02  1.647133e-02  5.581363e-03  5.054801e-03
##  [761] -6.557703e-03 -7.767662e-03  1.363973e-03  1.427769e-02 -2.418081e-03
##  [766]  1.819832e-02  3.270027e-03  2.301270e-02 -2.253797e-03 -2.991125e-02
##  [771]  1.207622e-02  1.873158e-03  9.966672e-03  7.709766e-03  2.598500e-03
##  [776]  3.280685e-02  1.625993e-03  1.439373e-02  5.219276e-03  1.623146e-02
##  [781]  3.861438e-02  2.123372e-02  3.600366e-03  1.337894e-02  3.187629e-02
##  [786]  2.744317e-02  5.016065e-03  6.945359e-02  4.004276e-02  6.432265e-02
##  [791] -8.231379e-02 -2.222791e-02 -1.879309e-01 -1.594491e-01  8.829088e-02
##  [796] -4.416791e-03  7.024146e-02  5.741264e-02 -1.005258e-03  5.063472e-02
##  [801]  8.855306e-03  3.417981e-02  2.042782e-02  2.099720e-02  3.775469e-02
##  [806]  5.556063e-03  3.736218e-03 -1.346561e-03  4.430428e-02  9.018224e-02
##  [811]  5.037513e-02  2.727601e-02 -3.614547e-02  5.310696e-02  4.517775e-02
##  [816]  4.077738e-02  3.272362e-02 -1.261612e-03 -4.656613e-02  3.182796e-02
##  [821] -4.027052e-03  9.564492e-02  9.936793e-02  3.623507e-02 -3.804907e-02
##  [826] -4.464457e-02  3.928935e-02  3.988421e-02  2.641992e-02  7.949359e-02
##  [831] -1.274707e-02  2.826361e-03  6.425468e-02  8.003002e-02  2.411099e-02
##  [836]  1.118276e-01  6.540735e-02 -6.841718e-03 -8.750231e-03 -8.783318e-04
##  [841]  1.264077e-03 -4.495753e-02 -1.235125e-01 -6.024069e-02 -1.465136e-02
##  [846]  3.218683e-03 -3.393316e-02  3.898582e-02 -7.331144e-03 -4.523821e-02
##  [851] -3.359462e-04  2.834409e-02 -1.248482e-01 -7.740810e-02  4.317631e-02
##  [856]  5.338504e-02 -1.069292e-02 -7.568585e-03  3.804925e-02  4.067247e-02
##  [861]  6.493314e-02 -1.211300e-02 -3.378259e-02 -2.081611e-02 -2.248239e-03
##  [866]  4.273058e-02  5.212959e-03 -4.984678e-02  4.568527e-02  2.208949e-02
##  [871] -1.665512e-02 -2.385893e-02 -1.786075e-03  1.613329e-03 -2.877193e-02
##  [876]  7.542516e-02  5.144290e-02  6.948621e-02  3.557614e-02 -1.577964e-02
##  [881]  1.015469e-02 -6.973152e-03 -5.335650e-02 -4.012849e-02 -6.516333e-02
##  [886] -2.409817e-02  1.311239e-02 -3.416647e-02 -6.051745e-02 -2.708777e-02
##  [891] -8.465658e-02  3.243353e-02 -3.489765e-03  6.594309e-03 -2.670611e-02
##  [896]  1.391400e-01  5.604973e-02 -3.491475e-02  4.335617e-02  2.132573e-02
##  [901]  2.914986e-02 -5.375900e-02 -3.671755e-02 -5.236379e-02 -3.749024e-03
##  [906] -1.249113e-01  4.686804e-02  3.698504e-02  5.709910e-02  1.623838e-02
##  [911] -9.346057e-02  4.316667e-02  8.094072e-04  4.421178e-02  1.276049e-02
##  [916]  1.003414e-02  5.853383e-02  7.601356e-02  2.479574e-02  1.502321e-02
##  [921] -1.494918e-02 -3.747079e-02  1.780053e-02  1.760261e-02 -5.747865e-02
##  [926] -6.508690e-02 -3.262799e-02 -7.983129e-02 -1.806984e-02  3.125303e-02
##  [931]  3.294768e-02  2.387051e-02  2.707049e-02 -2.916718e-03 -4.731836e-03
##  [936] -7.646974e-03  9.161886e-04 -5.335405e-02 -4.270356e-02 -5.418982e-04
##  [941]  4.868400e-02  1.821998e-02 -6.290235e-03  2.719435e-02  8.127192e-03
##  [946]  1.789163e-03 -2.470476e-02  1.902020e-02  4.639823e-03 -5.737484e-02
##  [951] -1.166978e-02  3.992005e-03 -5.651888e-03  4.767083e-04  8.724882e-03
##  [956] -5.284738e-02 -4.417280e-02  2.088911e-02  3.316477e-02 -2.822579e-04
##  [961]  1.276952e-02  2.278922e-02  1.345157e-02 -3.434570e-02 -2.367282e-02
##  [966]  1.662174e-02  1.153960e-03  1.932835e-02 -1.951854e-02 -5.145978e-02
##  [971] -4.882285e-02 -4.433700e-02 -2.293135e-02  1.304767e-02 -2.808784e-02
##  [976]  5.575945e-04 -1.435472e-02 -4.259757e-02 -5.816466e-02  2.255434e-02
##  [981] -3.409270e-02 -7.048140e-02 -3.120556e-02  5.840315e-03  3.794686e-02
##  [986]  1.567264e-02  7.980667e-03  3.867486e-02 -2.018860e-03  7.168722e-02
##  [991]  2.449481e-02 -4.152508e-02 -3.696586e-02 -6.872235e-02 -1.026006e-03
##  [996] -7.750335e-03 -7.885913e-03  4.330081e-02 -4.186811e-02 -2.087859e-02
## 
## $changePoints
## [1] 149 412 585 780 837
## 
## $detectionTimes
## [1] 253 477 662 803 980

The results at the monthly frequency are compared with those derived from the BSADF test (Figure 2 in the article).

#RENIXX
set.seed(004)
series <- ts(renixx_m[, c(2, 4)], frequency = 12, start = c(2005, 1))
results_r <- radf(series[, 1])
autoplot(results_r) + labs(title = "") + theme(
  plot.title = element_text(size = 20, face = "bold"),
  #Title size
  axis.title = element_text(size = 16),
  #Axis label size
  axis.text = element_text(size = 14),
  #Axis tick size
  legend.text = element_text(size = 12),
  #Legend text size
  legend.title = element_text(size = 14) #Legend title size
) + xlab('Time') + ylab('Value')
## Using `radf_crit` for `cv`.

d_adf_r <- datestamp(results_r)
## Using `radf_crit` for `cv`.
d_adf_r
## 
## ── Datestamp (min_duration = 0) ───────────────────────────────── Monte Carlo ──
## 
## series1 :
##        Start       Peak        End Duration   Signal Ongoing
## 1 2007-06-01 2007-07-01 2007-08-01        2 positive   FALSE
## 2 2007-10-01 2007-12-01 2008-01-01        3 positive   FALSE
## 3 2015-04-01 2015-04-01 2015-06-01        2 positive   FALSE
## 4 2019-12-01 2020-02-01 2020-03-01        3 positive   FALSE
## 5 2020-07-01 2021-01-01 2021-05-01       10 positive   FALSE
## 6 2021-11-01 2021-11-01 2021-12-01        1 positive   FALSE

As reported by the above results, there is a clear evidence that the BSADF is particularly effective in detecting explosive behaviors, a specific aspect of bubble dynamics. This highlights the need for a more comprehensive tool that encompasses the entire life cycle of a bubble, including its initiation, burst, and conclusion phases. This justify the adoption of new methodologies as the used KS-CPM.

Propagation phase

To analyze green bubbles, various variables will be taken into account. Two main groups stand out: economic/financial variables and Google Trends data.

#NOTE: This plot is for reference only and not included in the article.
#Brent Crude front month price
Oil <-
  read.csv("Economic variables/Oil_WTI.csv")
Oil$Date <- dmy(Oil$Exchange.Date)
Oil <- Oil[order(Oil$Date), ]
Oil <- Oil[, c(14, 2)]
rownames(Oil) <- NULL
Month <- as.Date(cut(Oil$Date, "month"))
Oil_m <- aggregate(Close ~ Month, Oil, mean)
#Transformation
Oil_m$Rate_O <- c(0, diff(log(Oil_m$Close)))
Oil_m$Rate_s <- scale(c(0, diff(Oil_m$Close)), scale = TRUE, center = TRUE)
rownames(Oil_m) <- NULL
ggplot(Oil_m, aes(Month, Close)) + geom_line() + xlab('Time') + ylab('Closed Value')

#NOTE: This plot is for reference only and not included in the article.
#MSCI_WORLD
MSCI <- read.csv("Economic variables/MSCI World.csv")
MSCI$Date <- dmy(MSCI$Exchange.Date)
MSCI <- MSCI[, c(9, 2)]
MSCI <- MSCI[order(MSCI$Date), ]
rownames(MSCI) <- NULL
Month <- as.Date(cut(MSCI$Date, "month"))
MSCI$Close <- as.numeric(gsub(",", "", MSCI$Close))
MSCI_m <- aggregate(Close ~ Month, MSCI, mean)
MSCI_m$Close <- as.numeric(MSCI_m$Close)
#Transformations
MSCI_m$Rate_ms <- c(0, diff(log(MSCI_m$Close)))
MSCI_m$Rate_s <-
  scale(c(0, diff(MSCI_m$Close)), scale = TRUE, center = TRUE)
ggplot(MSCI_m, aes(Month, Close)) + geom_line() + xlab('Time') + ylab('Closed Value')

The following section highlights the two bubbles identified by the KS-CPM in the renewable energy market. These findings are further supported by the behavior of the oil and MSCI ratio. The below chunk outlines the plot reported in Figure 3 in the article.

#Ratio
options(xts.message.period.apply.mean = FALSE)
renixx_m$Oil_ratio <- renixx_m$close / Oil_m$Close
renixx_m$MSCI_ratio <- renixx_m$close / MSCI_m$Close
#Final plot for green bubble detection (monthly)
renixx_m$Date <- as.POSIXct(renixx_m$Month, format = "%Y-%m-%d")
df_rect <-
  data.frame(
    xmin = c(as.POSIXct(d_adf_r$series1$Start)),
    xmax = c(as.POSIXct(d_adf_r$series1$End)),
    ymin = 0,
    ymax = Inf,
    Explosive = c("BSADF")
  )
g <- ggplot(renixx_m, aes(x = Date, y = close)) +
  geom_line(aes(color = "RENIXX")) +
  geom_line(aes(y = MSCI_ratio * 300, color = "Brent oil ratio"), linetype = "twodash") +
  geom_line(aes(y = Oil_ratio * 15, color = "MSCI ratio"), linetype = "twodash") +
  theme_light() +
  ylab('Value (€)') +
  xlab('Time') +
  ggtitle('Empirical Analysis (Monthly Time Series)') +
  theme(plot.title = element_text(hjust = 0.5)) +
  #Use annotate() instead of geom_segment() to avoid length warnings
  annotate(
    "segment",
    x = renixx_m$Date[36],
    xend = renixx_m$Date[36],
    y = 0,
    yend = renixx_m$close[36],
    linetype = "dashed",
    color = "red"
  ) +
  annotate(
    "segment",
    x = renixx_m$Date[95],
    xend = renixx_m$Date[95],
    y = 0,
    yend = renixx_m$close[95],
    linetype = "dashed",
    color = "red"
  ) +
  annotate(
    "segment",
    x = renixx_m$Date[124],
    xend = renixx_m$Date[124],
    y = 0,
    yend = renixx_m$close[124],
    linetype = "dashed",
    color = "red"
  ) +
  annotate(
    "segment",
    x = renixx_m$Date[178],
    xend = renixx_m$Date[178],
    y = 0,
    yend = renixx_m$close[178],
    linetype = "dashed",
    color = "red"
  ) +
  annotate(
    "segment",
    x = renixx_m$Date[193],
    xend = renixx_m$Date[193],
    y = 0,
    yend = renixx_m$close[193],
    linetype = "dashed",
    color = "red"
  ) +
  #Rectangle layer (fine as-is)
  geom_rect(
    data = df_rect,
    aes(
      xmin = xmin,
      ymin = ymin,
      xmax = xmax,
      ymax = ymax,
      fill = Explosive
    ),
    alpha = 0.2,
    inherit.aes = FALSE
  ) +
  scale_fill_manual(name = 'Test', values = c("azure4")) +
  #Highlight points
  geom_point(
    data = renixx_m[c(36, 95, 124, 178, 193), ],
    aes(x = Date, y = close),
    colour = "red",
    size = 2
  ) +
  #Manual legend and styling
  scale_color_manual(
    name = 'Time Series',
    breaks = c('RENIXX', 'Brent oil ratio', 'MSCI ratio'),
    values = c(
      'RENIXX' = 'darkgreen',
      'Brent oil ratio' = 'darkred',
      'MSCI ratio' = 'steelblue'
    )
  ) +
  theme(
    plot.title = element_text(hjust = 0, size = 20, face = "bold"),
    axis.title = element_text(size = 16),
    axis.text = element_text(size = 14),
    legend.text = element_text(size = 12),
    legend.title = element_text(size = 14)
  )
#Display the plot
grid.arrange(g)

The KS-CPM results are further validated through a robustness analysis using weekly frequency data, which confirms the previous findings based on low-frequency data (Figure A.1 in the article).

renixx_w$Date <- as.POSIXct(renixx_w$date, format = "%Y-%m-%d")
#Final plot for green bubble detection (weekly)
g <- ggplot(renixx_w, aes(x = Date, y = close)) +
  geom_line(aes(color = "RENIXX")) +
  theme_light() +
  ylab('Value (€)') +
  xlab('Time') +
  ggtitle('Empirical Analysis (Weekly Time Series)') +
  theme(plot.title = element_text(hjust = 0.5)) +
  #Use annotate() for vertical dashed lines
  annotate(
    "segment",
    x = renixx_w$Date[149],
    xend = renixx_w$Date[149],
    y = 0,
    yend = renixx_w$close[149],
    linetype = "dashed",
    color = "red"
  ) +
  annotate(
    "segment",
    x = renixx_w$Date[412],
    xend = renixx_w$Date[412],
    y = 0,
    yend = renixx_w$close[412],
    linetype = "dashed",
    color = "red"
  ) +
  annotate(
    "segment",
    x = renixx_w$Date[585],
    xend = renixx_w$Date[585],
    y = 0,
    yend = renixx_w$close[585],
    linetype = "dashed",
    color = "red"
  ) +
  annotate(
    "segment",
    x = renixx_w$Date[780],
    xend = renixx_w$Date[780],
    y = 0,
    yend = renixx_w$close[780],
    linetype = "dashed",
    color = "red"
  ) +
  annotate(
    "segment",
    x = renixx_w$Date[837],
    xend = renixx_w$Date[837],
    y = 0,
    yend = renixx_w$close[837],
    linetype = "dashed",
    color = "red"
  ) +
  #Highlight the points
  geom_point(
    data = renixx_w[c(149, 412, 585, 780, 837), ],
    aes(x = Date, y = close),
    colour = "red",
    size = 2
  ) +
  #Color legend and style
  scale_color_manual(
    name = 'Time Series',
    breaks = c('RENIXX', 'Brent oil ratio', 'MSCI ratio'),
    values = c(
      'RENIXX' = 'darkgreen',
      'Brent oil ratio' = 'darkred',
      'MSCI ratio' = 'steelblue'
    )
  ) +
  theme(
    plot.title = element_text(hjust = 0, size = 20, face = "bold"),
    axis.title = element_text(size = 16),
    axis.text = element_text(size = 14),
    legend.text = element_text(size = 12),
    legend.title = element_text(size = 14)
  )
grid.arrange(g)

Below, the GPR index is outlined.

#NOTE: This plot is for reference only and not included in the article.
#Geopolitical index
data_gpr_export <-
  read_excel("Economic variables/GPI.xls")
geo_i <- data_gpr_export %>% dplyr::select(6, 3)
colnames(geo_i)[1] <- "Date"
geo_i2 <-
  geo_i[geo_i$Date < "2024-03-01" &
          geo_i$Date >= "2005-01-01", ] #2005/01/01-2024/03/01
Month <- as.Date(cut(geo_i2$Date, "month"))
geo_m <- aggregate(geo_i2$GPRD ~ Month, geo_i2, mean)
colnames(geo_m)[2] <- "GPRH"
geo_m$GPRH <- as.numeric(geo_m$GPRH)
#Transformations
geo_m$Rate_g <- c(0, diff(log(geo_m$GPRH)))
geo_m$St <- scale(c(0, diff(geo_m$GPRH)), scale = TRUE, center = TRUE)
ggplot(geo_m, aes(Month, GPRH)) + geom_line() + xlab('Time') + ylab('Index Value')

Below, the FSI as well as the EPU index are outlined.

#NOTE: These plots are for reference only and not included in the article
#EPU
GEPUCURRENT <- read.csv("Economic variables/GEPUCURRENT.csv")
GEPUCURRENT$DATE <- ymd(GEPUCURRENT$DATE)
GEPUCURRENT <- GEPUCURRENT[97:326, ] #2005/01/01-2024/02/01
rownames(GEPUCURRENT) <- NULL
#Transformation
GEPUCURRENT$Rate <- c(0, diff(log(GEPUCURRENT$GEPUCURRENT)))
GEPUCURRENT$Rate_s <-
  scale(c(0, diff(GEPUCURRENT$GEPUCURRENT)), scale = TRUE, center = TRUE)
#OFR
Fs <- read.csv("Economic variables/FSI.csv")
Month <- as.Date(cut(ymd(Fs$Date), "month"))
Fs_m <- aggregate(OFR.FSI ~ Month, Fs, mean)
Fs_m <- Fs_m[61:290, ] #2005/01/01-2024/02/01
#Transformation
Fs_m$Rate <- c(0, diff(Fs_m$OFR.FSI))
Fs_m$Rate_s <-
  scale(c(0, diff(Fs_m$OFR.FSI)), scale = TRUE, center = TRUE)
ggplot(Fs_m, aes(Month, OFR.FSI)) + geom_line() + xlab('Time') + ylab('Index Value')

ggplot(GEPUCURRENT, aes(DATE, GEPUCURRENT)) + geom_line() + xlab('Time') + ylab('Index Value')

Finally, both economic and financial variables are grouped into a single dataset.

#NOTE: This plot is for reference only and not included in the article
#Financial/economic variables
Finance <- renixx_m[, c(1, 3)]
colnames(Finance)[2] <- "Renixx"
Finance$Oil_p <- Oil_m$Rate_O
Finance$MSCI <- MSCI_m$Rate_ms
Finance$EPU <- GEPUCURRENT$Rate
Finance$OFR <- Fs_m$Rate
Finance$Geo_i <- geo_m$Rate_g
#Financial dataset
Finance_2 <- renixx_m[, c(1, 2)]
colnames(Finance_2)[2] <- "Renixx"
Finance_2$Renixx <- scale(Finance_2$Renixx, center = TRUE, scale = TRUE)
Finance_2$Oil_p <- scale(Oil_m$Close, center = TRUE, scale = TRUE)
Finance_2$Geo_i <- scale(geo_m$GPRH, center = TRUE, scale = TRUE)
Finance_2$MSCI <- scale(MSCI_m$Close, center = TRUE, scale = TRUE)
Finance_2$EPU <- scale(GEPUCURRENT$GEPUCURRENT, center = TRUE, scale = TRUE)
Finance_2$OFR <- scale(Fs_m$OFR.FSI, center = TRUE, scale = TRUE)
#Plot
ggplot() +
  geom_line(data = Finance_2, aes(x = Month, y = Renixx, color = "Rennix")) +
  geom_line(data = Finance_2, aes(x = Month, y = MSCI, color = "Global MSCI")) +
  geom_line(data = Finance_2, aes(x = Month, y = Geo_i, color = "Geopolitical index")) +
  geom_line(data = Finance_2, aes(x = Month, y = Oil_p, color = "Oil price")) +
  geom_line(data = Finance_2, aes(x = Month, y = EPU, color = "EPU")) +
  geom_line(data = Finance_2, aes(x = Month, y = OFR, color = "OFR")) +
  xlab('Time') + ggtitle('Financial time series') +
  ylab('Monthly series') + scale_color_manual(
    name = 'Series',
    breaks = c(
      'Rennix',
      'Global MSCI',
      'Geopolitical index',
      'Oil price',
      'EPU',
      'OFR'
    ),
    values = c(
      'Rennix' = 'green',
      'Global MSCI' = 'red',
      'Geopolitical index' = 'blue',
      'Oil price' = 'brown',
      'EPU' = 'orange',
      'OFR' = 'grey'
    )
  )

Similarly, the text-based variables used in the article have been uploaded.

#NOTE: This plot is for reference only and not included in the article
#Green Energy
green_e <-
  read.csv("Google Trends/Green_energy.csv",
           header = TRUE,
           comment.char = "#")[13:242, ]
green_e$Time <- ymd(green_e$Time)
green_e$Rate <- c(0, diff(log(green_e$Absolute.Google.Search.Volume)))
Climate <- as.data.frame(green_e[, 1])
colnames(Climate)[1] <- "Date"
Climate$Green_e <- green_e[, 4]
#New technology
Tech_d <-
  read.csv("Google Trends/New_Technology.csv",
           header = TRUE,
           comment.char = "#")[13:242, ]
Tech_d$Rate <- c(0, diff(log(Tech_d$Absolute.Google.Search.Volume)))
Climate$Tech <- Tech_d[, 4]
#Energy index
Energy_i <-
  read.csv("Google Trends/Energy_index.csv",
           header = TRUE,
           comment.char = "#")[13:242, ]
Energy_i <- Energy_i[, c(1, 2, 3)]
Energy_i$Rate <-
  c(0, diff(log(Energy_i$Absolute.Google.Search.Volume)))
Climate$Energy_i <- Energy_i[, 4]
#Global warming
global_w <-
  read.csv("Google Trends/Global_warming.csv",
           header = TRUE,
           comment.char = "#")[13:242, ]
global_w$Rate <-
  c(0, diff(log(global_w$Absolute.Google.Search.Volume)))
Climate$Warming <- global_w[, 4]
#Natural disaster
Natural <-
  read.csv("Google Trends/Natural_disasters.csv",
           header = TRUE,
           comment.char = "#")[13:242, ]
Natural$Rate <- c(0, diff(log(Natural$Absolute.Google.Search.Volume)))
Climate$Natural <- Natural[, 4]
#Carbon price
Carbon_p <-
  read.csv("Google Trends/Carbon_price.csv",
           header = TRUE,
           comment.char = "#")[13:242, ]
Carbon_p$Rate <-
  c(0, diff(log(Carbon_p$Absolute.Google.Search.Volume)))
Climate$Carbon_p <- Carbon_p[, 4]
#Carbon tax
Tax <-
  read.csv("Google Trends/Carbon_tax.csv",
           header = TRUE,
           comment.char = "#")[13:242, ]
Tax$Rate <- c(0, diff(log(Tax$Absolute.Google.Search.Volume)))
Climate$Tax <- Tax[, 4]
#Energy shares
Energy_s <-
  read.csv("Google Trends/Energy_shares.csv",
           header = TRUE,
           comment.char = "#")[13:242, ]
Energy_s$Rate <-
  c(0, diff(log(Energy_s$Absolute.Google.Search.Volume)))
Climate$Energy_s <- Energy_s[, 4]
#Text dataset
Climate_2 <- as.data.frame(green_e[, 1])
colnames(Climate_2)[1] <- "Date"
#Scaling
Climate_2$Green_e <-
  scale(as.numeric(green_e[, 3]), center = TRUE, scale = TRUE)
Climate_2$Tech <-
  scale(as.numeric(Tech_d[, 3]), center = TRUE, scale = TRUE)
Climate_2$Warming <-
  scale(as.numeric(global_w[, 3]), center = TRUE, scale = TRUE)
Climate_2$Energy_i <-
  scale(as.numeric(Energy_i[, 3]), center = TRUE, scale = TRUE)
Climate_2$Natural <-
  scale(as.numeric(Natural[, 3]), center = TRUE, scale = TRUE)
Climate_2$Carbon_p <-
  scale(as.numeric(Carbon_p[, 3]), center = TRUE, scale = TRUE)
Climate_2$Energy_s <-
  scale(as.numeric(Energy_s[, 3]), center = TRUE, scale = TRUE)
Climate_2$Energy_T <-
  scale(as.numeric(Energy_s[, 3] + Energy_i[, 3]) / 2,
        center = TRUE,
        scale = TRUE)
Climate_2$Energy_s <-
  scale(as.numeric(Energy_s[, 3]), center = TRUE, scale = TRUE)
Climate_2$Tax <- scale(as.numeric(Tax[, 3]), center = TRUE, scale = TRUE)
#Plot
ggplot() +
  geom_line(data = Climate_2, aes(x = Date, y = Green_e, color = "Green Energy")) +
  geom_line(data = Climate_2, aes(x = Date, y = Energy_i, color = "Energy Index")) +
  geom_line(data = Climate_2, aes(x = Date, y = Warming, color = "Global Warming")) +
  geom_line(data = Climate_2, aes(x = Date, y = Natural, color = "Natural Disasters")) +
  geom_line(data = Climate_2, aes(x = Date, y = Carbon_p, color = "Carbon price")) +
  geom_line(data = Climate_2, aes(x = Date, y = Energy_s, color = "Energy shares")) +
  geom_line(data = Climate_2, aes(x = Date, y = Tax, color = "Carbon Tax")) +
  geom_line(data = Climate_2, aes(x = Date, y = Tech, color = "New Technology")) +
  xlab('Time') + ggtitle('Text-time series') +
  ylab('Monthly series') + scale_color_manual(
    name = 'Series',
    breaks = c(
      'Green Energy',
      'Global Warming',
      'Natural Disasters',
      'Carbon price',
      "New Technology",
      'Energy Index',
      'Green Bond',
      'Carbon Emissions',
      'Energy shares',
      'Carbon Tax'
    ),
    values = c(
      'Green Energy' = 'green',
      'Global Warming' = 'red',
      'Natural Disasters' = 'gold',
      'Carbon price' = 'brown',
      'New Technology' = 'purple',
      'Energy Index' = 'blue',
      'Carbon Emissions' = 'orange',
      'Energy shares' = 'violet',
      'Carbon Tax' = 'grey'
    )
  )

The final dataset is built as below.

#Final dataset (log-diff)
Finance$Merge <- format(as.Date(Finance$Month), "%Y-%m")
Climate$Merge <- format(as.Date(Climate$Date), "%Y-%m")
Data <- merge(Finance, Climate, by = "Merge")
Data <- Data[, -c(1, 9)]
rownames(Data) <- NULL
#Final dataset (scaled)
Climate_2 <-
  Climate_2[c(
    'Date',
    'Green_e',
    'Tech',
    'Warming',
    'Natural',
    'Carbon_p',
    'Energy_i',
    'Energy_s',
    'Tax'
  )]
Finance_2$Merge <- format(as.Date(Finance_2$Month), "%Y-%m")
Climate_2$Merge <- format(as.Date(Climate_2$Date), "%Y-%m")
#Origianl dataset
Data_or <- merge(Finance_2, Climate_2, by = "Merge")
Data_or <- Data_or[, -c(1, 9)]
rownames(Data_or) <- NULL

Before proceeding, both the Augmented Dickey-Fuller (ADF) test are conducted to detect the presence of a unit root.

#NOTE: These were included in a previous version of the article and are retained for reference.
#ADF for text data
summary(ur.df(
  Climate[, "Carbon_p"],
  #lag=12 for seasonality
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41381 -0.07046 -0.01482  0.05959  0.47782 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.008949   0.007835   1.142  0.25462    
## z.lag.1     -1.674943   0.108639 -15.417  < 2e-16 ***
## z.diff.lag   0.205206   0.065164   3.149  0.00186 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.118 on 225 degrees of freedom
## Multiple R-squared:  0.7083, Adjusted R-squared:  0.7057 
## F-statistic: 273.1 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -15.4175 118.8527 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
summary(ur.df(
  Climate[, "Green_e"],
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34956 -0.07668  0.00021  0.06728  0.38783 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.007327   0.007507   0.976     0.33    
## z.lag.1     -1.568111   0.099976 -15.685   <2e-16 ***
## z.diff.lag   0.278004   0.063799   4.358    2e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1131 on 225 degrees of freedom
## Multiple R-squared:  0.6452, Adjusted R-squared:  0.6421 
## F-statistic: 204.6 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -15.6848 123.0077 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
summary(ur.df(
  Climate[, "Tech"],
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.225374 -0.067981 -0.006719  0.057915  0.265361 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.00543    0.00642  -0.846    0.399    
## z.lag.1     -1.24245    0.08901 -13.959  < 2e-16 ***
## z.diff.lag   0.28840    0.06410   4.500 1.09e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09673 on 225 degrees of freedom
## Multiple R-squared:  0.5241, Adjusted R-squared:  0.5199 
## F-statistic: 123.9 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -13.9588 97.4254 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
summary(ur.df(
  Climate[, "Warming"],
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4295 -0.1267  0.0110  0.1194  0.7071 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.00585    0.01242  -0.471   0.6382    
## z.lag.1     -1.02124    0.08733 -11.694   <2e-16 ***
## z.diff.lag   0.15827    0.06577   2.407   0.0169 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1875 on 225 degrees of freedom
## Multiple R-squared:  0.4553, Adjusted R-squared:  0.4505 
## F-statistic: 94.04 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -11.6939 68.3743 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
summary(ur.df(
  Climate[, "Natural"],
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.61849 -0.15003  0.01427  0.13392  1.09063 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.006564   0.016519  -0.397    0.691    
## z.lag.1     -1.374322   0.087720 -15.667  < 2e-16 ***
## z.diff.lag   0.360640   0.061775   5.838 1.83e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2493 on 225 degrees of freedom
## Multiple R-squared:  0.5736, Adjusted R-squared:  0.5699 
## F-statistic: 151.4 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -15.6672 122.746 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
summary(ur.df(
  Climate[, "Tax"],
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.11273 -0.21206 -0.00117  0.21038  1.23944 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.009043   0.024992   0.362 0.717805    
## z.lag.1     -1.489789   0.099825 -14.924  < 2e-16 ***
## z.diff.lag   0.245664   0.064678   3.798 0.000188 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3772 on 225 degrees of freedom
## Multiple R-squared:  0.6228, Adjusted R-squared:  0.6195 
## F-statistic: 185.8 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -14.924 111.368 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
summary(ur.df(
  Climate[, "Energy_s"],
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.24993 -0.11629 -0.00469  0.12850  0.82393 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.00909    0.01428   0.636 0.525160    
## z.lag.1     -1.79177    0.10992 -16.301  < 2e-16 ***
## z.diff.lag   0.23776    0.06459   3.681 0.000291 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2155 on 225 degrees of freedom
## Multiple R-squared:  0.7405, Adjusted R-squared:  0.7382 
## F-statistic:   321 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -16.3005 132.8539 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
summary(ur.df(
  Climate[, "Energy_i"],
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.45124 -0.09009 -0.00872  0.08769  0.51330 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.002114   0.009444   0.224  0.82305    
## z.lag.1     -1.411003   0.100610 -14.024  < 2e-16 ***
## z.diff.lag   0.186673   0.065297   2.859  0.00465 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1426 on 225 degrees of freedom
## Multiple R-squared:  0.6102, Adjusted R-squared:  0.6068 
## F-statistic: 176.1 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -14.0245 98.3442 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
#ADF for economic data
summary(ur.df(
  Finance[, "Oil_p"],
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47455 -0.04535  0.00878  0.04822  0.23000 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.002470   0.005456   0.453   0.6512    
## z.lag.1     -0.737127   0.076217  -9.671   <2e-16 ***
## z.diff.lag   0.114039   0.066237   1.722   0.0865 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0823 on 225 degrees of freedom
## Multiple R-squared:  0.3394, Adjusted R-squared:  0.3336 
## F-statistic: 57.81 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -9.6714 46.768 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
summary(ur.df(
  Finance[, "Renixx"],
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.43442 -0.04634  0.00149  0.04643  0.22351 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.002994   0.005575   0.537   0.5918    
## z.lag.1     -0.787127   0.078432 -10.036   <2e-16 ***
## z.diff.lag   0.122234   0.066193   1.847   0.0661 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08401 on 225 degrees of freedom
## Multiple R-squared:  0.3603, Adjusted R-squared:  0.3546 
## F-statistic: 63.37 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -10.0357 50.3618 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
summary(ur.df(
  Finance[, "MSCI"],
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.243032 -0.013601  0.005433  0.022454  0.108286 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.005151   0.002566   2.007   0.0459 *  
## z.lag.1     -0.958193   0.090121 -10.632   <2e-16 ***
## z.diff.lag   0.052153   0.066709   0.782   0.4352    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03808 on 225 degrees of freedom
## Multiple R-squared:  0.4559, Adjusted R-squared:  0.451 
## F-statistic: 94.25 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -10.6323 56.5239 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
summary(ur.df(
  Finance[, "EPU"],
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.43837 -0.11260 -0.02077  0.10251  0.60925 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.00843    0.01157   0.729   0.4669    
## z.lag.1     -1.44543    0.10392 -13.910   <2e-16 ***
## z.diff.lag   0.16543    0.06595   2.509   0.0128 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1744 on 225 degrees of freedom
## Multiple R-squared:  0.6296, Adjusted R-squared:  0.6263 
## F-statistic: 191.3 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -13.9096 96.7406 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
summary(ur.df(
  Finance[, "OFR"],
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3733 -0.5323 -0.1483  0.3311 11.7428 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.009295   0.098257   0.095   0.9247    
## z.lag.1     -0.911928   0.082130 -11.103   <2e-16 ***
## z.diff.lag   0.167023   0.065697   2.542   0.0117 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.484 on 225 degrees of freedom
## Multiple R-squared:  0.4079, Adjusted R-squared:  0.4027 
## F-statistic: 77.52 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -11.1034 61.643 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
summary(ur.df(
  Finance[, "Geo_i"],
  type = "drift",
  lags = 1,
  selectlags = "AIC"
))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49065 -0.11058 -0.01727  0.09786  0.67808 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.00205    0.01236   0.166   0.8684    
## z.lag.1     -1.46444    0.10522 -13.918   <2e-16 ***
## z.diff.lag   0.14898    0.06595   2.259   0.0248 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1866 on 225 degrees of freedom
## Multiple R-squared:  0.6455, Adjusted R-squared:  0.6423 
## F-statistic: 204.8 on 2 and 225 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -13.918 96.8558 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
#Arch test due to no log-transformation
archtest(Finance[, "OFR"])
## 
##  Engle's LM ARCH Test
## 
## data:  Finance[, "OFR"]
## statistic = 1.7834, lag = 1, p-value = 0.1817
## alternative hypothesis: ARCH effects of order 1 are present

The propagation effects of both the Renixx index and oil prices are evaluated using econometric techniques. Specifically, a cointegration analysis is conducted on the variables under consideration (e.g. Johansen).

#Function to obtain log-likelihood from VECM
logLik_cajorls <- function(cajorl_model) {
  #Extract residuals
  resids <- residuals(cajorl_model$rlm)
  #Sample size (T) and number of equations (K)
  T <- nrow(resids)
  K <- ncol(resids)
  #Residual covariance matrix
  Sigma <- crossprod(resids) / T
  #Log determinant of Sigma
  log_det_Sigma <- log(det(Sigma))
  #Log-likelihood formula for multivariate normal residuals
  logLik <- -T / 2 * (K * log(2 * pi) + log_det_Sigma + K)
  return(logLik)
}

All series are integrated of order one, I(1), and the Johansen procedure indicates one cointegrating relationship (r=1) at the 1% significance level and two at the 5% significance level (Table 3 in the article).

#VAR (12 time frequency of the series)
all_series <-
  ts(
    Data_or[, c('Renixx', 'OFR', 'Oil_p')],
    start = c(2005, 1),
    end = c(2024, 2),
    frequency = 12
  )
#Johansen procedure
jotest = ca.jo(
  all_series,
  type = "trace",
  K = 2,
  ecdet = "none",
  spec = "longrun"
)
summary(jotest)
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: trace statistic , with linear trend 
## 
## Eigenvalues (lambda):
## [1] 0.10210818 0.05986997 0.02200688
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 2 |  5.07  6.50  8.18 11.65
## r <= 1 | 19.15 15.66 17.95 23.52
## r = 0  | 43.71 28.71 31.52 37.22
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##            Renixx.l2    OFR.l2   Oil_p.l2
## Renixx.l2  1.0000000 1.0000000  1.0000000
## OFR.l2    -1.7967895 0.1675129  0.4801047
## Oil_p.l2  -0.2248293 5.0244366 -0.3037410
## 
## Weights W:
## (This is the loading matrix)
## 
##           Renixx.l2       OFR.l2    Oil_p.l2
## Renixx.d 0.00179232 -0.006868744 -0.01603928
## OFR.d    0.04144494  0.011258129 -0.00101930
## Oil_p.d  0.01390801 -0.009160432  0.01940998

The VECM results indicate that both Renixx.dl1 and OFR.dl1 significantly affect OFR, suggesting that these variables are important short-term drivers of the dependent variable (Table 4 in the article).

vecm_model <- cajorls(jotest, r = 2)
summary(vecm_model$rlm)[2]
## Response OFR.d :
## 
## Call:
## lm(formula = OFR.d ~ ect1 + ect2 + constant + Renixx.dl1 + OFR.dl1 + 
##     Oil_p.dl1 - 1, data = data.mat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.74983 -0.12657 -0.04008  0.06820  2.79772 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## ect1        0.052703   0.012540   4.203 3.82e-05 ***
## ect2       -0.072582   0.021254  -3.415 0.000758 ***
## constant    0.001395   0.021326   0.065 0.947914    
## Renixx.dl1  0.249939   0.125811   1.987 0.048194 *  
## OFR.dl1     0.191546   0.077994   2.456 0.014821 *  
## Oil_p.dl1  -0.028551   0.085874  -0.332 0.739839    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3215 on 222 degrees of freedom
## Multiple R-squared:  0.1493, Adjusted R-squared:  0.1263 
## F-statistic: 6.492 on 6 and 222 DF,  p-value: 2.516e-06
cat("Log-likelihood:", logLik_cajorls(vecm_model), "\n")
## Log-likelihood: 31.74526

The following chunks are dedicated to analyzing the dynamics of bubbles in the Renixx index.

Each time series is divided according to the three stages detected by the KS-CPM.

#Codification
breack
## [1] "2007-12-01" "2012-11-01" "2015-04-01" "2019-10-01" "2021-01-01"
Data_2 <- renixx_m[c(2:230), c('Month', 'close')] #For lagged values
colnames(Data_2)[2] <- "Renixx"
#Dummy codification
Data_2$dummy  <- cut(
  Data_2$Month,
  breaks = as.Date(c(
    "2005-01-01", "2012-11-01", "2019-10-01", '2024-03-01'
  )),
  labels = c('Clean-Tech', 'No-bubble', 'Climate')
)
#Dataset lagged values of text variables only (scaled values)
Data_2$dummy <- relevel(Data_2$dummy, ref = 'No-bubble')
Data_2$Energy_ilag <- lag(Data_or$Energy_i)[2:230]
Data_2$Warminglag <- lag(Data_or$Warming)[2:230]
Data_2$Naturallag <- lag(Data_or$Natural)[2:230]
Data_2$Green_elag <- lag(Data_or$Green_e)[2:230]
Data_2$Carbon_plag <- lag(Data_or$Carbon_p)[2:230]
Data_2$Techlag <- lag(Data_or$Tech)[2:230]
Data_2$Energy_slag <- lag(Data_or$Energy_s)[2:230]
Data_2$Taxlag <- lag(Data_or$Tax)[2:230]
rownames(Data_2) <- NULL
#Info
str(Data_2)
## 'data.frame':    229 obs. of  11 variables:
##  $ Month      : Date, format: "2005-02-01" "2005-03-01" ...
##  $ Renixx     : num  421 451 445 444 503 ...
##  $ dummy      : Factor w/ 3 levels "No-bubble","Clean-Tech",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Energy_ilag: num  0.1 0.821 1.723 0.581 0.461 ...
##  $ Warminglag : num  0.565 0.801 0.506 0.388 0.506 ...
##  $ Naturallag : num  7.45 3.39 2.11 2.01 1.71 ...
##  $ Green_elag : num  -1.76 -1.32 -1.47 -1.18 -1.47 ...
##  $ Carbon_plag: num  -1.19 -1.3 -1.09 -1.25 -1.14 ...
##  $ Techlag    : num  3.61 3.69 3.44 3.2 2.79 ...
##  $ Energy_slag: num  -1.075 -0.526 -2.051 -0.709 -1.197 ...
##  $ Taxlag     : num  -0.933 -0.737 -0.835 -0.835 -0.442 ...

The plot below illustrates the distribution of text variables based on the three different stages identified by the KS-CPM (Figure 4 in the article).

#Density plot for the bubble's stages
p1 <-
  ggplot(Data_2,
         aes(
           x = Green_elag,
           group = dummy,
           color = dummy,
           fill = dummy
         )) + labs(fill = 'Stage') +
  guides(color = "none") + geom_density(alpha = 0.4) + ggtitle('Green Energy') + ylab('Density') + xlab('Standardized Value') +
  xlim(c(-2, 4)) + ylim(c(0, 1.5))
p2 <-
  ggplot(Data_2,
         aes(
           x = Naturallag,
           group = dummy,
           color = dummy,
           fill = dummy
         )) + geom_density(alpha = 0.4) +
  guides(color = "none") + geom_density(alpha = 0.4) + ggtitle('Natural Disasters') + labs(fill =
                                                                                             'Stage') + ylab('Density') + xlab('Standardized Value') +
  xlim(c(-2, 4)) + ylim(c(0, 1.5))
p3 <-
  ggplot(Data_2,
         aes(
           x = Carbon_plag,
           group = dummy,
           color = dummy,
           fill = dummy
         )) + geom_density(alpha = 0.4) + ggtitle('Carbon Price') +
  guides(color = "none") + geom_density(alpha = 0.4) + labs(fill = 'Stage') + xlab('Standardized Value') + ylab('Density') +
  xlim(c(-2, 4)) + ylim(c(0, 1.5))
p4 <-
  ggplot(Data_2,
         aes(
           x = Energy_ilag,
           group = dummy,
           color = dummy,
           fill = dummy
         )) + geom_density(alpha = 0.4) +
  guides(color = "none") + geom_density(alpha = 0.4) + xlab('Standardized Value') + ggtitle('Energy Index') + labs(fill =
                                                                                                                     'Stage') + ylab('Density') +
  xlim(c(-2, 4)) + ylim(c(0, 1.8))
p5 <-
  ggplot(Data_2,
         aes(
           x = Warminglag,
           group = dummy,
           color = dummy,
           fill = dummy
         )) + geom_density(alpha = 0.4) +
  guides(color = "none") + geom_density(alpha = 0.4) + xlab('Standardized Value') + ggtitle('Global Warming') + labs(fill =
                                                                                                                       'Stage') + ylab('Density') +
  xlim(c(-2, 4)) + ylim(c(0, 2.8))
p6 <-
  ggplot(Data_2, aes(
    x = Techlag,
    group = dummy,
    color = dummy,
    fill = dummy
  )) + geom_density(alpha = 0.4) +
  guides(color = "none") + geom_density(alpha = 0.4) + xlab('Standardized Value') +
  ggtitle('New Technology') + labs(fill =
                                     'Stage') + ylab('Density') +
  xlim(c(-2, 4)) + ylim(c(0, 1.5))
p7 <-
  ggplot(Data_2,
         aes(
           x = Energy_slag,
           group = dummy,
           color = dummy,
           fill = dummy
         )) + geom_density(alpha = 0.4) +
  guides(color = "none") + geom_density(alpha = 0.4) + xlab('Standardized Value') + ggtitle('Energy Shares') + labs(fill =
                                                                                                                      'Stage') + ylab('Density') +
  xlim(c(-2, 4)) + ylim(c(0, 1.5))
p8 <-
  ggplot(Data_2, aes(
    x = Taxlag,
    group = dummy,
    color = dummy,
    fill = dummy
  )) + geom_density(alpha = 0.4) +
  guides(color = "none") + geom_density(alpha = 0.4) + xlab('Standardized Value') + ggtitle('Carbon Tax') + labs(fill =
                                                                                                                   'Stage') + ylab('Density') +
  xlim(c(-2, 4)) + ylim(c(0, 1.5))
suppressWarnings({
  grid.arrange(
    p1,
    p2,
    p3,
    p4,
    p5,
    p6,
    p7,
    p8,
    ncol = 4,
    top = textGrob("RENIXX Index", gp = gpar(fontsize = 20, font = 3))
  )
})

Particularly, the subplot related to the term green energy delineates how the phase of the Climate bubble in the Renixx index (depicted in blue) shares similarities with that of the Clean-Tech bubble (illustrated in green), differing from the No-bubble phase (represented in red). This indicates a comparable level of research interest during speculative phases for the term green energy. These findings are consistent with existing literature, reinforcing the central narrative of the Clean-Tech bubble as one of . Moreover, the terms and exhibit higher research volumes compared to previous levels (e.g. Climate bubble phase).

This is further supported by the results of the non-parametric KS test applied to the three previously mentioned Google Trends and their respective distributions (Table 2 in the article).

#KSTEST Energy Index
x <- Data_2[Data_2$dummy == 'Climate', 'Energy_ilag']
y <- Data_2[Data_2$dummy == 'No-bubble', 'Energy_ilag']
z <- Data_2[Data_2$dummy == 'Clean-Tech', 'Energy_ilag']
cat("\033[34mKS test (p-value) for Energy Index (No-bubble):",
    "\033[39m\n")
## [34mKS test (p-value) for Energy Index (No-bubble): [39m
ks.test(x, y, alternative = 'l')
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D^- = 0.9434, p-value < 2.2e-16
## alternative hypothesis: the CDF of x lies below that of y
cat("\033[34mKS test (p-value) for Energy Index (Clean-Tech):",
    "\033[39m\n")
## [34mKS test (p-value) for Energy Index (Clean-Tech): [39m
ks.test(x, z, alternative = 'l')
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  x and z
## D^- = 0.60945, p-value = 4.359e-13
## alternative hypothesis: the CDF of x lies below that of y
#KSTEST Energy Shares
x <- Data_2[Data_2$dummy == 'Climate', 'Energy_slag']
y <- Data_2[Data_2$dummy == 'No-bubble', 'Energy_slag']
z <- Data_2[Data_2$dummy == 'Clean-Tech', 'Energy_slag']
cat("\033[34mKS test (p-value) for Energy Shares (No-bubble):",
    "\033[39m\n")
## [34mKS test (p-value) for Energy Shares (No-bubble): [39m
ks.test(x, y, alternative = 'l')
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D^- = 0.96908, p-value < 2.2e-16
## alternative hypothesis: the CDF of x lies below that of y
cat("\033[34mKS test (p-value) for Energy Shares (Clean-Tech):",
    "\033[39m\n")
## [34mKS test (p-value) for Energy Shares (Clean-Tech): [39m
ks.test(x, z, alternative = 'l')
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  x and z
## D^- = 0.85474, p-value < 2.2e-16
## alternative hypothesis: the CDF of x lies below that of y
#KSTEST Green energy
x <- Data_2[Data_2$dummy == 'Climate', 'Green_elag']
y <- Data_2[Data_2$dummy == 'No-bubble', 'Green_elag']
z <- Data_2[Data_2$dummy == 'Clean-Tech', 'Green_elag']
cat("\033[34mKS test (p-value) for Green Energy (No-bubble):",
    "\033[39m\n")
## [34mKS test (p-value) for Green Energy (No-bubble): [39m
ks.test(z, y, alternative = 'l')
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  z and y
## D^- = 0.56341, p-value = 1.971e-14
## alternative hypothesis: the CDF of x lies below that of y
cat("\033[34mKS test (p-value) for Green Energy (Climate):",
    "\033[39m\n")
## [34mKS test (p-value) for Green Energy (Climate): [39m
ks.test(z, x, alternative = 'l')
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  z and x
## D^- = -1.1796e-16, p-value = 1
## alternative hypothesis: the CDF of x lies below that of y

Forecasting Phase

Considering the limited sample size given by monthly time seires data, in the forecasting phase other than Renixx other 3 indexes are taken into account for the forecasting phase.

#iShares_GCE
iShares_Clean_Energy <- read.csv("Economic variables/iShares_Clean_Energy.csv")
iShares_Clean_Energy$Date <- dmy(iShares_Clean_Energy$Exchange.Date)
iShares_Clean_Energy <- iShares_Clean_Energy[, c(17, 2)]
iShares_Clean_Energy <- iShares_Clean_Energy[order(iShares_Clean_Energy$Date), ]
rownames(iShares_Clean_Energy) <- NULL
iShares_Clean_Energy$Close <- as.numeric(gsub(",", "", iShares_Clean_Energy$Close))
Month <- as.Date(cut(iShares_Clean_Energy$Date, "month"))
iShares_Clean_Energy_m <- aggregate(Close ~ Month, iShares_Clean_Energy, mean)
iShares_Clean_Energy_m$Close <- as.numeric(iShares_Clean_Energy_m$Close)
#SP_GCE
SP_Clean_Energy <- read.csv("Economic variables/S&P_Clean_Energy.csv")
SP_Clean_Energy$Date <- dmy(SP_Clean_Energy$Exchange.Date)
SP_Clean_Energy <- SP_Clean_Energy[, c(5, 2)]
SP_Clean_Energy <- SP_Clean_Energy[order(SP_Clean_Energy$Date), ]
rownames(SP_Clean_Energy) <- NULL
Month <- as.Date(cut(SP_Clean_Energy$Date, "month"))
SP_Clean_Energy$Close <- as.numeric(gsub(",", "", SP_Clean_Energy$Close))
SP_Clean_Energy_m <- aggregate(Close ~ Month, SP_Clean_Energy, mean)
SP_Clean_Energy_m$Close <- as.numeric(SP_Clean_Energy_m$Close)
#MSCI_GAE
MSCI_Alternative_Energy_Price <- read.csv("Economic variables/MSCI_Alternative_Energy_Price.csv")
MSCI_Alternative_Energy_Price$Date <- dmy(MSCI_Alternative_Energy_Price$Exchange.Date)
MSCI_Alternative_Energy_Price <- MSCI_Alternative_Energy_Price[, c(9, 2)]
MSCI_Alternative_Energy_Price <- MSCI_Alternative_Energy_Price[order(MSCI_Alternative_Energy_Price$Date), ]
rownames(MSCI_Alternative_Energy_Price) <- NULL
Month <- as.Date(cut(MSCI_Alternative_Energy_Price$Date, "month"))
MSCI_Alternative_Energy_Price$Close <- as.numeric(gsub(",", "", MSCI_Alternative_Energy_Price$Close))
MSCI_Alternative_Energy_Price_m <- aggregate(Close ~ Month, MSCI_Alternative_Energy_Price, mean)
MSCI_Alternative_Energy_Price_m$Close <- as.numeric(MSCI_Alternative_Energy_Price_m$Close)

Because of the escalating utilization of Google throughout the years, the chosen time frame extends from 2015 to nowadays. This decision aims to capture a heightened correlation between text variables and the Renixx index, as depicted in the above plot. In the one-step ahead exercise, four models are taken into account.

Below are the summary statistics for the considered dataset (Table 5 in the article).

#Dataset
data_pred <-
  data.frame(renixx_m$close, green_e[, 1])[121:230, ] #Date
colnames(data_pred) <-
  c('Renixx', 'Date')
#iShares_GCE
data_pred$iShares_GCE <- iShares_Clean_Energy_m$Close[91:200]
#SP_GCE
data_pred$SP_GCE <- SP_Clean_Energy_m$Close[73:182]
#MSCI_GAE
data_pred$MSCI_GAE <- MSCI_Alternative_Energy_Price_m$Close[52:161]
options(digits = 9) #For 3-digits results
summary(data_pred[, c(1, 5, 3, 4)])
##      Renixx            MSCI_GAE         iShares_GCE           SP_GCE        
##  Min.   : 387.821   Min.   : 43.7995   Min.   : 4.59050   Min.   : 597.911  
##  1st Qu.: 452.972   1st Qu.: 49.1302   1st Qu.: 5.32000   1st Qu.: 701.396  
##  Median : 560.122   Median : 54.5514   Median : 6.30478   Median : 869.168  
##  Mean   : 900.370   Mean   : 63.8702   Mean   : 7.97172   Mean   :1119.184  
##  3rd Qu.:1475.525   3rd Qu.: 80.4228   3rd Qu.:10.88143   3rd Qu.:1590.921  
##  Max.   :2169.222   Max.   :119.3667   Max.   :18.04600   Max.   :2574.594
#Dataset
data_pred <-
  data.frame(
    log(renixx_m$close),
    log(MSCI_m$Close),
    log(geo_m$GPRH),
    log(GEPUCURRENT$GEPUCURRENT),
    log(Oil_m$Close),
    log(Energy_i[, 3]),
    log(Energy_s[, 3]),
    log(green_e[, 3]),
    log(global_w[, 3]),
    log(Natural[, 3]),
    log(Carbon_p[, 3]),
    log(Tax[, 3]),
    log(Tech_d[, 3]),
    green_e[, 1] #Date
  )[121:230, ] #2015/01/01 2024/02/01 and log transformation
colnames(data_pred) <-
  c(
    'Renixx',
    'MSCI',
    'Geop',
    'EPU',
    'Oil_p',
    'Energy_i',
    'Energy_s',
    'Green_e',
    'Global_w',
    'Natural_d',
    'Carbon_p',
    'Carbon_t',
    'Tech',
    'Date'
  )
rownames(data_pred) <- NULL
#iShares_GCE
data_pred$iShares_GCE <- log(iShares_Clean_Energy_m$Close[91:200])
#SP_GCE
data_pred$SP_GCE <- log(SP_Clean_Energy_m$Close[73:182])
#MSCI_GAE
data_pred$MSCI_GAE <- log(MSCI_Alternative_Energy_Price_m$Close[52:161])
#Iterative short-term forecasts
iteration <- 37 #Bubble burst
models <- 4
sample <- length(data_pred[, 1])
#RENIXX
predictions_renixx <- as.data.frame(matrix(nrow = iteration, ncol = models))
colnames(predictions_renixx) <-
  c('ARIMA', 'BSTS', 'ARIMAX_t1', 'ARIMAX_t')
up_uni_renixx <- as.data.frame(matrix(nrow = iteration, ncol = models))
colnames(up_uni_renixx) <- c('ARIMA', 'BSTS', 'ARIMAX_t1', 'ARIMAX_t')
lower_uni_renixx <- as.data.frame(matrix(nrow = iteration, ncol = models))
colnames(lower_uni_renixx) <- c(c('ARIMA', 'BSTS', 'ARIMAX_t1', 'ARIMAX_t'))
#MSCI_GAE
predictions_MSCI_GAE <- as.data.frame(matrix(nrow = iteration, ncol = models))
colnames(predictions_MSCI_GAE) <-
  c('ARIMA', 'BSTS', 'ARIMAX_t1', 'ARIMAX_t')
up_uni_MSCI_GAE <- as.data.frame(matrix(nrow = iteration, ncol = models))
colnames(up_uni_MSCI_GAE) <- c('ARIMA', 'BSTS', 'ARIMAX_t1', 'ARIMAX_t')
lower_uni_MSCI_GAE <- as.data.frame(matrix(nrow = iteration, ncol = models))
colnames(lower_uni_MSCI_GAE) <- c(c('ARIMA', 'BSTS', 'ARIMAX_t1', 'ARIMAX_t'))
predictions_SP_GCE <- as.data.frame(matrix(nrow = iteration, ncol = models))
#SP_GCE
colnames(predictions_SP_GCE) <-
  c('ARIMA', 'BSTS', 'ARIMAX_t1', 'ARIMAX_t')
predictions_iShares_GCE <- as.data.frame(matrix(nrow = iteration, ncol = models))
up_uni_SP_GCE <- as.data.frame(matrix(nrow = iteration, ncol = models))
colnames(up_uni_SP_GCE) <- c('ARIMA', 'BSTS', 'ARIMAX_t1', 'ARIMAX_t')
lower_uni_SP_GCE <- as.data.frame(matrix(nrow = iteration, ncol = models))
colnames(lower_uni_SP_GCE) <- c(c('ARIMA', 'BSTS', 'ARIMAX_t1', 'ARIMAX_t'))
#iShares_GCE
colnames(predictions_iShares_GCE) <-
  c('ARIMA', 'BSTS', 'ARIMAX_t1', 'ARIMAX_t')
up_uni_iShares_GCE <- as.data.frame(matrix(nrow = iteration, ncol = models))
colnames(up_uni_iShares_GCE) <- c('ARIMA', 'BSTS', 'ARIMAX_t1', 'ARIMAX_t')
lower_uni_iShares_GCE <- as.data.frame(matrix(nrow = iteration, ncol = models))
colnames(lower_uni_iShares_GCE) <- c(c('ARIMA', 'BSTS', 'ARIMAX_t1', 'ARIMAX_t'))

Consistent with existing literature, the only exogenous variable considered for inclusion in the ARIMAX model is the Google Trends energy index. This is further supported by a multivariate Granger causality test with lagged values (p=1), where the Google Trends energy index is the only variable showing a significance level below 1% across all four considered indexes. However, the EPU index could also be a potential candidate for inclusion in the model (Table 7 in the article).

Note: refer to the fourth column of the first table in each multivariate test.

#Granger causality (multivariate). In the paper, the reported value is the p-value corresponding to the F-statistic (first table - fourth  column).
cat("\033[34mRENIXX:", "\033[39m\n")
## [34mRENIXX: [39m
vm = VAR(data_pred[1:110, c(1, 3:4, 6:13)], p = 1, type = 'both')
granger_causality(vm)
## 
## Granger Causality Test (Multivariate)
## 
## F test and Wald χ² test based on VAR(1) model:
## ────────────────────────────────────────────────────────────────────
##                               F df1 df2     p     Chisq df     p    
## ────────────────────────────────────────────────────────────────────
##  ---------------------                                              
##  Renixx <= Geop            0.00   1  96  .973      0.00  1  .973    
##  Renixx <= EPU             3.70   1  96  .057 .    3.70  1  .054 .  
##  Renixx <= Energy_i       14.51   1  96 <.001 *** 14.51  1 <.001 ***
##  Renixx <= Energy_s        1.56   1  96  .214      1.56  1  .211    
##  Renixx <= Green_e         0.21   1  96  .651      0.21  1  .650    
##  Renixx <= Global_w        0.50   1  96  .479      0.50  1  .478    
##  Renixx <= Natural_d       0.03   1  96  .874      0.03  1  .874    
##  Renixx <= Carbon_p        0.05   1  96  .818      0.05  1  .818    
##  Renixx <= Carbon_t        0.21   1  96  .648      0.21  1  .647    
##  Renixx <= Tech            2.97   1  96  .088 .    2.97  1  .085 .  
##  Renixx <= ALL             3.03  10  96  .002 **  30.27 10 <.001 ***
##  ---------------------                                              
##  Geop <= Renixx            1.42   1  96  .236      1.42  1  .233    
##  Geop <= EPU               1.45   1  96  .232      1.45  1  .229    
##  Geop <= Energy_i          1.24   1  96  .269      1.24  1  .266    
##  Geop <= Energy_s          1.23   1  96  .270      1.23  1  .268    
##  Geop <= Green_e           0.37   1  96  .542      0.37  1  .541    
##  Geop <= Global_w          0.09   1  96  .767      0.09  1  .766    
##  Geop <= Natural_d         1.16   1  96  .283      1.16  1  .281    
##  Geop <= Carbon_p          1.26   1  96  .265      1.26  1  .262    
##  Geop <= Carbon_t          0.09   1  96  .771      0.09  1  .770    
##  Geop <= Tech              2.51   1  96  .116      2.51  1  .113    
##  Geop <= ALL               1.53  10  96  .141     15.30 10  .122    
##  ---------------------                                              
##  EPU <= Renixx             0.41   1  96  .525      0.41  1  .523    
##  EPU <= Geop               0.43   1  96  .513      0.43  1  .511    
##  EPU <= Energy_i           0.22   1  96  .643      0.22  1  .641    
##  EPU <= Energy_s           1.23   1  96  .271      1.23  1  .268    
##  EPU <= Green_e            2.63   1  96  .108      2.63  1  .105    
##  EPU <= Global_w           1.14   1  96  .288      1.14  1  .286    
##  EPU <= Natural_d          4.16   1  96  .044 *    4.16  1  .041 *  
##  EPU <= Carbon_p           1.19   1  96  .277      1.19  1  .275    
##  EPU <= Carbon_t           0.54   1  96  .464      0.54  1  .462    
##  EPU <= Tech               3.02   1  96  .086 .    3.02  1  .082 .  
##  EPU <= ALL                2.03  10  96  .038 *   20.34 10  .026 *  
##  ---------------------                                              
##  Energy_i <= Renixx        3.25   1  96  .075 .    3.25  1  .072 .  
##  Energy_i <= Geop          0.22   1  96  .639      0.22  1  .638    
##  Energy_i <= EPU           0.56   1  96  .457      0.56  1  .455    
##  Energy_i <= Energy_s     11.06   1  96  .001 **  11.06  1 <.001 ***
##  Energy_i <= Green_e       0.46   1  96  .501      0.46  1  .499    
##  Energy_i <= Global_w      0.23   1  96  .633      0.23  1  .632    
##  Energy_i <= Natural_d     0.12   1  96  .727      0.12  1  .726    
##  Energy_i <= Carbon_p      0.29   1  96  .592      0.29  1  .591    
##  Energy_i <= Carbon_t      1.46   1  96  .230      1.46  1  .227    
##  Energy_i <= Tech          2.29   1  96  .134      2.29  1  .130    
##  Energy_i <= ALL           2.20  10  96  .024 *   22.04 10  .015 *  
##  ---------------------                                              
##  Energy_s <= Renixx        0.11   1  96  .739      0.11  1  .738    
##  Energy_s <= Geop          0.37   1  96  .544      0.37  1  .543    
##  Energy_s <= EPU           0.00   1  96  .988      0.00  1  .988    
##  Energy_s <= Energy_i      0.68   1  96  .410      0.68  1  .408    
##  Energy_s <= Green_e       0.24   1  96  .627      0.24  1  .626    
##  Energy_s <= Global_w      0.54   1  96  .465      0.54  1  .464    
##  Energy_s <= Natural_d     0.24   1  96  .625      0.24  1  .624    
##  Energy_s <= Carbon_p      3.48   1  96  .065 .    3.48  1  .062 .  
##  Energy_s <= Carbon_t      0.31   1  96  .577      0.31  1  .576    
##  Energy_s <= Tech          4.00   1  96  .048 *    4.00  1  .045 *  
##  Energy_s <= ALL           1.49  10  96  .156     14.87 10  .137    
##  ---------------------                                              
##  Green_e <= Renixx         0.00   1  96  .953      0.00  1  .953    
##  Green_e <= Geop           1.42   1  96  .237      1.42  1  .234    
##  Green_e <= EPU            0.07   1  96  .791      0.07  1  .790    
##  Green_e <= Energy_i       0.20   1  96  .657      0.20  1  .656    
##  Green_e <= Energy_s       1.12   1  96  .293      1.12  1  .291    
##  Green_e <= Global_w       5.30   1  96  .023 *    5.30  1  .021 *  
##  Green_e <= Natural_d      0.57   1  96  .453      0.57  1  .451    
##  Green_e <= Carbon_p       2.20   1  96  .142      2.20  1  .138    
##  Green_e <= Carbon_t       0.05   1  96  .818      0.05  1  .817    
##  Green_e <= Tech           1.66   1  96  .200      1.66  1  .197    
##  Green_e <= ALL            1.76  10  96  .079 .   17.60 10  .062 .  
##  ---------------------                                              
##  Global_w <= Renixx        1.97   1  96  .164      1.97  1  .161    
##  Global_w <= Geop          1.28   1  96  .261      1.28  1  .258    
##  Global_w <= EPU           0.03   1  96  .873      0.03  1  .872    
##  Global_w <= Energy_i      4.36   1  96  .039 *    4.36  1  .037 *  
##  Global_w <= Energy_s      1.62   1  96  .206      1.62  1  .203    
##  Global_w <= Green_e       0.92   1  96  .339      0.92  1  .337    
##  Global_w <= Natural_d     0.01   1  96  .931      0.01  1  .931    
##  Global_w <= Carbon_p      0.02   1  96  .885      0.02  1  .885    
##  Global_w <= Carbon_t      2.30   1  96  .133      2.30  1  .129    
##  Global_w <= Tech          6.61   1  96  .012 *    6.61  1  .010 *  
##  Global_w <= ALL           1.68  10  96  .096 .   16.84 10  .078 .  
##  ---------------------                                              
##  Natural_d <= Renixx       0.03   1  96  .872      0.03  1  .872    
##  Natural_d <= Geop         1.28   1  96  .261      1.28  1  .258    
##  Natural_d <= EPU          0.17   1  96  .681      0.17  1  .680    
##  Natural_d <= Energy_i     3.05   1  96  .084 .    3.05  1  .081 .  
##  Natural_d <= Energy_s     2.17   1  96  .144      2.17  1  .141    
##  Natural_d <= Green_e      2.20   1  96  .141      2.20  1  .138    
##  Natural_d <= Global_w     0.00   1  96  .951      0.00  1  .951    
##  Natural_d <= Carbon_p     0.07   1  96  .795      0.07  1  .794    
##  Natural_d <= Carbon_t     0.99   1  96  .323      0.99  1  .321    
##  Natural_d <= Tech         5.40   1  96  .022 *    5.40  1  .020 *  
##  Natural_d <= ALL          1.62  10  96  .112     16.23 10  .093 .  
##  ---------------------                                              
##  Carbon_p <= Renixx        2.92   1  96  .091 .    2.92  1  .088 .  
##  Carbon_p <= Geop          0.95   1  96  .332      0.95  1  .329    
##  Carbon_p <= EPU           1.67   1  96  .199      1.67  1  .196    
##  Carbon_p <= Energy_i      0.52   1  96  .474      0.52  1  .473    
##  Carbon_p <= Energy_s      0.22   1  96  .640      0.22  1  .639    
##  Carbon_p <= Green_e       0.67   1  96  .415      0.67  1  .413    
##  Carbon_p <= Global_w      0.36   1  96  .550      0.36  1  .548    
##  Carbon_p <= Natural_d     0.46   1  96  .498      0.46  1  .497    
##  Carbon_p <= Carbon_t      2.57   1  96  .112      2.57  1  .109    
##  Carbon_p <= Tech          4.47   1  96  .037 *    4.47  1  .034 *  
##  Carbon_p <= ALL           1.62  10  96  .112     16.21 10  .094 .  
##  ---------------------                                              
##  Carbon_t <= Renixx        1.13   1  96  .290      1.13  1  .287    
##  Carbon_t <= Geop          0.50   1  96  .483      0.50  1  .481    
##  Carbon_t <= EPU           0.32   1  96  .570      0.32  1  .569    
##  Carbon_t <= Energy_i      3.33   1  96  .071 .    3.33  1  .068 .  
##  Carbon_t <= Energy_s      3.41   1  96  .068 .    3.41  1  .065 .  
##  Carbon_t <= Green_e       0.41   1  96  .522      0.41  1  .521    
##  Carbon_t <= Global_w      6.59   1  96  .012 *    6.59  1  .010 *  
##  Carbon_t <= Natural_d     0.00   1  96  .969      0.00  1  .969    
##  Carbon_t <= Carbon_p      0.21   1  96  .647      0.21  1  .646    
##  Carbon_t <= Tech          2.85   1  96  .095 .    2.85  1  .091 .  
##  Carbon_t <= ALL           2.81  10  96  .004 **  28.11 10  .002 ** 
##  ---------------------                                              
##  Tech <= Renixx            0.45   1  96  .503      0.45  1  .502    
##  Tech <= Geop              8.16   1  96  .005 **   8.16  1  .004 ** 
##  Tech <= EPU               0.37   1  96  .546      0.37  1  .544    
##  Tech <= Energy_i          1.47   1  96  .228      1.47  1  .225    
##  Tech <= Energy_s          0.21   1  96  .648      0.21  1  .647    
##  Tech <= Green_e           2.00   1  96  .161      2.00  1  .157    
##  Tech <= Global_w          1.22   1  96  .271      1.22  1  .268    
##  Tech <= Natural_d         0.64   1  96  .427      0.64  1  .425    
##  Tech <= Carbon_p          1.87   1  96  .174      1.87  1  .171    
##  Tech <= Carbon_t          0.15   1  96  .696      0.15  1  .695    
##  Tech <= ALL               2.69  10  96  .006 **  26.92 10  .003 ** 
## ────────────────────────────────────────────────────────────────────
cat("\033[34mMSCI_GAE :", "\033[39m\n")
## [34mMSCI_GAE : [39m
vm = VAR(data_pred[1:110, c(17, 3:4, 6:13)], p = 1, type = 'both')
granger_causality(vm)
## 
## Granger Causality Test (Multivariate)
## 
## F test and Wald χ² test based on VAR(1) model:
## ────────────────────────────────────────────────────────────────────
##                               F df1 df2     p     Chisq df     p    
## ────────────────────────────────────────────────────────────────────
##  ---------------------                                              
##  MSCI_GAE <= Geop          0.56   1  96  .456      0.56  1  .454    
##  MSCI_GAE <= EPU           6.75   1  96  .011 *    6.75  1  .009 ** 
##  MSCI_GAE <= Energy_i      9.46   1  96  .003 **   9.46  1  .002 ** 
##  MSCI_GAE <= Energy_s      0.97   1  96  .327      0.97  1  .325    
##  MSCI_GAE <= Green_e       0.08   1  96  .777      0.08  1  .776    
##  MSCI_GAE <= Global_w      0.25   1  96  .620      0.25  1  .619    
##  MSCI_GAE <= Natural_d     0.00   1  96  .961      0.00  1  .961    
##  MSCI_GAE <= Carbon_p      0.11   1  96  .737      0.11  1  .736    
##  MSCI_GAE <= Carbon_t      0.05   1  96  .825      0.05  1  .825    
##  MSCI_GAE <= Tech          2.68   1  96  .105      2.68  1  .101    
##  MSCI_GAE <= ALL           2.99  10  96  .003 **  29.88 10 <.001 ***
##  ---------------------                                              
##  Geop <= MSCI_GAE          0.79   1  96  .378      0.79  1  .375    
##  Geop <= EPU               1.24   1  96  .268      1.24  1  .265    
##  Geop <= Energy_i          0.88   1  96  .350      0.88  1  .348    
##  Geop <= Energy_s          1.40   1  96  .240      1.40  1  .237    
##  Geop <= Green_e           0.52   1  96  .474      0.52  1  .473    
##  Geop <= Global_w          0.07   1  96  .794      0.07  1  .794    
##  Geop <= Natural_d         1.04   1  96  .310      1.04  1  .307    
##  Geop <= Carbon_p          0.61   1  96  .435      0.61  1  .433    
##  Geop <= Carbon_t          0.04   1  96  .837      0.04  1  .836    
##  Geop <= Tech              2.63   1  96  .108      2.63  1  .105    
##  Geop <= ALL               1.46  10  96  .168     14.57 10  .149    
##  ---------------------                                              
##  EPU <= MSCI_GAE           1.33   1  96  .252      1.33  1  .249    
##  EPU <= Geop               0.57   1  96  .451      0.57  1  .449    
##  EPU <= Energy_i           0.49   1  96  .485      0.49  1  .483    
##  EPU <= Energy_s           1.14   1  96  .289      1.14  1  .286    
##  EPU <= Green_e            2.51   1  96  .116      2.51  1  .113    
##  EPU <= Global_w           1.48   1  96  .226      1.48  1  .223    
##  EPU <= Natural_d          3.49   1  96  .065 .    3.49  1  .062 .  
##  EPU <= Carbon_p           0.80   1  96  .374      0.80  1  .372    
##  EPU <= Carbon_t           0.24   1  96  .622      0.24  1  .621    
##  EPU <= Tech               3.20   1  96  .077 .    3.20  1  .074 .  
##  EPU <= ALL                2.14  10  96  .028 *   21.45 10  .018 *  
##  ---------------------                                              
##  Energy_i <= MSCI_GAE      2.64   1  96  .107      2.64  1  .104    
##  Energy_i <= Geop          0.18   1  96  .674      0.18  1  .673    
##  Energy_i <= EPU           0.47   1  96  .495      0.47  1  .493    
##  Energy_i <= Energy_s     11.83   1  96 <.001 *** 11.83  1 <.001 ***
##  Energy_i <= Green_e       0.33   1  96  .569      0.33  1  .568    
##  Energy_i <= Global_w      0.13   1  96  .717      0.13  1  .716    
##  Energy_i <= Natural_d     0.13   1  96  .721      0.13  1  .720    
##  Energy_i <= Carbon_p      1.29   1  96  .259      1.29  1  .256    
##  Energy_i <= Carbon_t      1.42   1  96  .237      1.42  1  .234    
##  Energy_i <= Tech          2.18   1  96  .143      2.18  1  .140    
##  Energy_i <= ALL           2.13  10  96  .029 *   21.32 10  .019 *  
##  ---------------------                                              
##  Energy_s <= MSCI_GAE      0.10   1  96  .751      0.10  1  .750    
##  Energy_s <= Geop          0.36   1  96  .548      0.36  1  .546    
##  Energy_s <= EPU           0.00   1  96  .994      0.00  1  .994    
##  Energy_s <= Energy_i      0.74   1  96  .391      0.74  1  .389    
##  Energy_s <= Green_e       0.23   1  96  .636      0.23  1  .635    
##  Energy_s <= Global_w      0.56   1  96  .457      0.56  1  .456    
##  Energy_s <= Natural_d     0.24   1  96  .623      0.24  1  .622    
##  Energy_s <= Carbon_p      4.15   1  96  .044 *    4.15  1  .042 *  
##  Energy_s <= Carbon_t      0.31   1  96  .579      0.31  1  .577    
##  Energy_s <= Tech          4.03   1  96  .048 *    4.03  1  .045 *  
##  Energy_s <= ALL           1.49  10  96  .157     14.85 10  .137    
##  ---------------------                                              
##  Green_e <= MSCI_GAE       0.51   1  96  .475      0.51  1  .473    
##  Green_e <= Geop           1.76   1  96  .188      1.76  1  .185    
##  Green_e <= EPU            0.19   1  96  .667      0.19  1  .666    
##  Green_e <= Energy_i       0.02   1  96  .890      0.02  1  .890    
##  Green_e <= Energy_s       1.15   1  96  .287      1.15  1  .284    
##  Green_e <= Global_w       5.86   1  96  .017 *    5.86  1  .015 *  
##  Green_e <= Natural_d      0.87   1  96  .352      0.87  1  .350    
##  Green_e <= Carbon_p       2.55   1  96  .114      2.55  1  .111    
##  Green_e <= Carbon_t       0.25   1  96  .618      0.25  1  .617    
##  Green_e <= Tech           1.84   1  96  .178      1.84  1  .175    
##  Green_e <= ALL            1.82  10  96  .067 .   18.20 10  .052 .  
##  ---------------------                                              
##  Global_w <= MSCI_GAE      4.54   1  96  .036 *    4.54  1  .033 *  
##  Global_w <= Geop          1.63   1  96  .205      1.63  1  .202    
##  Global_w <= EPU           0.11   1  96  .744      0.11  1  .744    
##  Global_w <= Energy_i      5.98   1  96  .016 *    5.98  1  .014 *  
##  Global_w <= Energy_s      1.94   1  96  .167      1.94  1  .164    
##  Global_w <= Green_e       1.08   1  96  .301      1.08  1  .299    
##  Global_w <= Natural_d     0.03   1  96  .874      0.03  1  .874    
##  Global_w <= Carbon_p      0.13   1  96  .721      0.13  1  .720    
##  Global_w <= Carbon_t      3.45   1  96  .066 .    3.45  1  .063 .  
##  Global_w <= Tech          7.08   1  96  .009 **   7.08  1  .008 ** 
##  Global_w <= ALL           1.98  10  96  .044 *   19.80 10  .031 *  
##  ---------------------                                              
##  Natural_d <= MSCI_GAE     0.44   1  96  .509      0.44  1  .508    
##  Natural_d <= Geop         1.52   1  96  .221      1.52  1  .218    
##  Natural_d <= EPU          0.08   1  96  .772      0.08  1  .772    
##  Natural_d <= Energy_i     4.07   1  96  .046 *    4.07  1  .044 *  
##  Natural_d <= Energy_s     2.23   1  96  .138      2.23  1  .135    
##  Natural_d <= Green_e      2.49   1  96  .118      2.49  1  .115    
##  Natural_d <= Global_w     0.00   1  96  .954      0.00  1  .954    
##  Natural_d <= Carbon_p     0.04   1  96  .852      0.04  1  .852    
##  Natural_d <= Carbon_t     0.59   1  96  .443      0.59  1  .442    
##  Natural_d <= Tech         5.64   1  96  .019 *    5.64  1  .018 *  
##  Natural_d <= ALL          1.67  10  96  .099 .   16.71 10  .081 .  
##  ---------------------                                              
##  Carbon_p <= MSCI_GAE      0.13   1  96  .717      0.13  1  .717    
##  Carbon_p <= Geop          1.53   1  96  .219      1.53  1  .216    
##  Carbon_p <= EPU           0.86   1  96  .355      0.86  1  .353    
##  Carbon_p <= Energy_i      0.00   1  96  .983      0.00  1  .983    
##  Carbon_p <= Energy_s      0.13   1  96  .716      0.13  1  .715    
##  Carbon_p <= Green_e       1.24   1  96  .269      1.24  1  .266    
##  Carbon_p <= Global_w      0.19   1  96  .661      0.19  1  .660    
##  Carbon_p <= Natural_d     0.95   1  96  .332      0.95  1  .329    
##  Carbon_p <= Carbon_t      1.16   1  96  .284      1.16  1  .281    
##  Carbon_p <= Tech          3.79   1  96  .054 .    3.79  1  .052 .  
##  Carbon_p <= ALL           1.31  10  96  .239     13.05 10  .221    
##  ---------------------                                              
##  Carbon_t <= MSCI_GAE      1.27   1  96  .263      1.27  1  .260    
##  Carbon_t <= Geop          0.50   1  96  .479      0.50  1  .477    
##  Carbon_t <= EPU           0.32   1  96  .572      0.32  1  .570    
##  Carbon_t <= Energy_i      3.49   1  96  .065 .    3.49  1  .062 .  
##  Carbon_t <= Energy_s      3.16   1  96  .079 .    3.16  1  .075 .  
##  Carbon_t <= Green_e       0.47   1  96  .495      0.47  1  .494    
##  Carbon_t <= Global_w      6.92   1  96  .010 **   6.92  1  .009 ** 
##  Carbon_t <= Natural_d     0.01   1  96  .927      0.01  1  .927    
##  Carbon_t <= Carbon_p      0.01   1  96  .903      0.01  1  .903    
##  Carbon_t <= Tech          2.87   1  96  .093 .    2.87  1  .090 .  
##  Carbon_t <= ALL           2.83  10  96  .004 **  28.28 10  .002 ** 
##  ---------------------                                              
##  Tech <= MSCI_GAE          1.25   1  96  .265      1.25  1  .263    
##  Tech <= Geop              8.74   1  96  .004 **   8.74  1  .003 ** 
##  Tech <= EPU               0.26   1  96  .614      0.26  1  .613    
##  Tech <= Energy_i          2.05   1  96  .156      2.05  1  .152    
##  Tech <= Energy_s          0.26   1  96  .611      0.26  1  .610    
##  Tech <= Green_e           2.18   1  96  .143      2.18  1  .139    
##  Tech <= Global_w          1.55   1  96  .217      1.55  1  .214    
##  Tech <= Natural_d         0.42   1  96  .519      0.42  1  .518    
##  Tech <= Carbon_p          2.79   1  96  .098 .    2.79  1  .095 .  
##  Tech <= Carbon_t          0.03   1  96  .858      0.03  1  .857    
##  Tech <= ALL               2.79  10  96  .004 **  27.94 10  .002 ** 
## ────────────────────────────────────────────────────────────────────
cat("\033[34miShares_GCE:", "\033[39m\n")
## [34miShares_GCE: [39m
vm = VAR(data_pred[1:110, c(15, 3:4, 6:13)], p = 1, type = 'both')
granger_causality(vm)
## 
## Granger Causality Test (Multivariate)
## 
## F test and Wald χ² test based on VAR(1) model:
## ───────────────────────────────────────────────────────────────────────
##                                  F df1 df2     p     Chisq df     p    
## ───────────────────────────────────────────────────────────────────────
##  ------------------------                                              
##  iShares_GCE <= Geop          0.06   1  96  .800      0.06  1  .800    
##  iShares_GCE <= EPU           7.40   1  96  .008 **   7.40  1  .007 ** 
##  iShares_GCE <= Energy_i     12.90   1  96 <.001 *** 12.90  1 <.001 ***
##  iShares_GCE <= Energy_s      2.63   1  96  .108      2.63  1  .105    
##  iShares_GCE <= Green_e       0.27   1  96  .605      0.27  1  .604    
##  iShares_GCE <= Global_w      0.73   1  96  .396      0.73  1  .394    
##  iShares_GCE <= Natural_d     0.53   1  96  .468      0.53  1  .466    
##  iShares_GCE <= Carbon_p      0.43   1  96  .513      0.43  1  .512    
##  iShares_GCE <= Carbon_t      0.26   1  96  .614      0.26  1  .613    
##  iShares_GCE <= Tech          1.19   1  96  .278      1.19  1  .276    
##  iShares_GCE <= ALL           3.38  10  96 <.001 *** 33.80 10 <.001 ***
##  ------------------------                                              
##  Geop <= iShares_GCE          2.07   1  96  .153      2.07  1  .150    
##  Geop <= EPU                  1.79   1  96  .183      1.79  1  .180    
##  Geop <= Energy_i             1.46   1  96  .230      1.46  1  .227    
##  Geop <= Energy_s             1.33   1  96  .251      1.33  1  .248    
##  Geop <= Green_e              0.18   1  96  .670      0.18  1  .669    
##  Geop <= Global_w             0.05   1  96  .819      0.05  1  .818    
##  Geop <= Natural_d            1.30   1  96  .257      1.30  1  .254    
##  Geop <= Carbon_p             1.13   1  96  .290      1.13  1  .288    
##  Geop <= Carbon_t             0.11   1  96  .737      0.11  1  .737    
##  Geop <= Tech                 2.15   1  96  .146      2.15  1  .143    
##  Geop <= ALL                  1.60  10  96  .117     16.04 10  .099 .  
##  ------------------------                                              
##  EPU <= iShares_GCE           1.35   1  96  .249      1.35  1  .246    
##  EPU <= Geop                  0.64   1  96  .425      0.64  1  .423    
##  EPU <= Energy_i              0.49   1  96  .485      0.49  1  .483    
##  EPU <= Energy_s              1.21   1  96  .274      1.21  1  .271    
##  EPU <= Green_e               1.93   1  96  .168      1.93  1  .165    
##  EPU <= Global_w              1.36   1  96  .246      1.36  1  .243    
##  EPU <= Natural_d             3.70   1  96  .057 .    3.70  1  .054 .  
##  EPU <= Carbon_p              1.29   1  96  .258      1.29  1  .255    
##  EPU <= Carbon_t              0.35   1  96  .555      0.35  1  .553    
##  EPU <= Tech                  3.44   1  96  .067 .    3.44  1  .064 .  
##  EPU <= ALL                   2.15  10  96  .028 *   21.47 10  .018 *  
##  ------------------------                                              
##  Energy_i <= iShares_GCE      3.36   1  96  .070 .    3.36  1  .067 .  
##  Energy_i <= Geop             0.29   1  96  .593      0.29  1  .592    
##  Energy_i <= EPU              0.74   1  96  .392      0.74  1  .390    
##  Energy_i <= Energy_s        11.57   1  96 <.001 *** 11.57  1 <.001 ***
##  Energy_i <= Green_e          0.68   1  96  .413      0.68  1  .411    
##  Energy_i <= Global_w         0.18   1  96  .670      0.18  1  .669    
##  Energy_i <= Natural_d        0.13   1  96  .723      0.13  1  .722    
##  Energy_i <= Carbon_p         0.54   1  96  .464      0.54  1  .462    
##  Energy_i <= Carbon_t         1.64   1  96  .203      1.64  1  .200    
##  Energy_i <= Tech             2.61   1  96  .109      2.61  1  .106    
##  Energy_i <= ALL              2.22  10  96  .023 *   22.17 10  .014 *  
##  ------------------------                                              
##  Energy_s <= iShares_GCE      0.43   1  96  .512      0.43  1  .511    
##  Energy_s <= Geop             0.48   1  96  .490      0.48  1  .488    
##  Energy_s <= EPU              0.02   1  96  .889      0.02  1  .889    
##  Energy_s <= Energy_i         0.49   1  96  .484      0.49  1  .483    
##  Energy_s <= Green_e          0.36   1  96  .549      0.36  1  .548    
##  Energy_s <= Global_w         0.62   1  96  .431      0.62  1  .429    
##  Energy_s <= Natural_d        0.33   1  96  .570      0.33  1  .569    
##  Energy_s <= Carbon_p         3.59   1  96  .061 .    3.59  1  .058 .  
##  Energy_s <= Carbon_t         0.45   1  96  .506      0.45  1  .505    
##  Energy_s <= Tech             3.67   1  96  .058 .    3.67  1  .055 .  
##  Energy_s <= ALL              1.52  10  96  .143     15.24 10  .124    
##  ------------------------                                              
##  Green_e <= iShares_GCE       1.10   1  96  .297      1.10  1  .294    
##  Green_e <= Geop              2.04   1  96  .157      2.04  1  .154    
##  Green_e <= EPU               0.34   1  96  .561      0.34  1  .559    
##  Green_e <= Energy_i          0.00   1  96  .987      0.00  1  .986    
##  Green_e <= Energy_s          1.10   1  96  .297      1.10  1  .295    
##  Green_e <= Global_w          6.07   1  96  .015 *    6.07  1  .014 *  
##  Green_e <= Natural_d         1.00   1  96  .319      1.00  1  .316    
##  Green_e <= Carbon_p          1.88   1  96  .173      1.88  1  .170    
##  Green_e <= Carbon_t          0.33   1  96  .566      0.33  1  .564    
##  Green_e <= Tech              2.11   1  96  .150      2.11  1  .147    
##  Green_e <= ALL               1.89  10  96  .056 .   18.90 10  .042 *  
##  ------------------------                                              
##  Global_w <= iShares_GCE      3.82   1  96  .054 .    3.82  1  .051 .  
##  Global_w <= Geop             1.72   1  96  .193      1.72  1  .190    
##  Global_w <= EPU              0.18   1  96  .673      0.18  1  .672    
##  Global_w <= Energy_i         5.54   1  96  .021 *    5.54  1  .019 *  
##  Global_w <= Energy_s         1.77   1  96  .186      1.77  1  .183    
##  Global_w <= Green_e          1.50   1  96  .223      1.50  1  .220    
##  Global_w <= Natural_d        0.00   1  96  .965      0.00  1  .965    
##  Global_w <= Carbon_p         0.01   1  96  .931      0.01  1  .931    
##  Global_w <= Carbon_t         2.84   1  96  .095 .    2.84  1  .092 .  
##  Global_w <= Tech             7.55   1  96  .007 **   7.55  1  .006 ** 
##  Global_w <= ALL              1.90  10  96  .055 .   18.97 10  .041 *  
##  ------------------------                                              
##  Natural_d <= iShares_GCE     0.51   1  96  .477      0.51  1  .476    
##  Natural_d <= Geop            1.59   1  96  .210      1.59  1  .207    
##  Natural_d <= EPU             0.05   1  96  .828      0.05  1  .827    
##  Natural_d <= Energy_i        4.18   1  96  .044 *    4.18  1  .041 *  
##  Natural_d <= Energy_s        2.17   1  96  .144      2.17  1  .140    
##  Natural_d <= Green_e         2.66   1  96  .106      2.66  1  .103    
##  Natural_d <= Global_w        0.00   1  96  .977      0.00  1  .977    
##  Natural_d <= Carbon_p        0.11   1  96  .735      0.11  1  .735    
##  Natural_d <= Carbon_t        0.68   1  96  .412      0.68  1  .410    
##  Natural_d <= Tech            5.81   1  96  .018 *    5.81  1  .016 *  
##  Natural_d <= ALL             1.68  10  96  .097 .   16.79 10  .079 .  
##  ------------------------                                              
##  Carbon_p <= iShares_GCE      1.95   1  96  .166      1.95  1  .163    
##  Carbon_p <= Geop             0.95   1  96  .331      0.95  1  .329    
##  Carbon_p <= EPU              1.63   1  96  .205      1.63  1  .202    
##  Carbon_p <= Energy_i         0.27   1  96  .607      0.27  1  .606    
##  Carbon_p <= Energy_s         0.16   1  96  .692      0.16  1  .691    
##  Carbon_p <= Green_e          0.55   1  96  .460      0.55  1  .458    
##  Carbon_p <= Global_w         0.35   1  96  .557      0.35  1  .555    
##  Carbon_p <= Natural_d        0.57   1  96  .453      0.57  1  .451    
##  Carbon_p <= Carbon_t         2.10   1  96  .151      2.10  1  .148    
##  Carbon_p <= Tech             4.58   1  96  .035 *    4.58  1  .032 *  
##  Carbon_p <= ALL              1.51  10  96  .147     15.11 10  .128    
##  ------------------------                                              
##  Carbon_t <= iShares_GCE      0.27   1  96  .604      0.27  1  .603    
##  Carbon_t <= Geop             0.37   1  96  .547      0.37  1  .546    
##  Carbon_t <= EPU              0.44   1  96  .510      0.44  1  .508    
##  Carbon_t <= Energy_i         4.65   1  96  .034 *    4.65  1  .031 *  
##  Carbon_t <= Energy_s         3.20   1  96  .077 .    3.20  1  .073 .  
##  Carbon_t <= Green_e          0.48   1  96  .490      0.48  1  .488    
##  Carbon_t <= Global_w         6.19   1  96  .015 *    6.19  1  .013 *  
##  Carbon_t <= Natural_d        0.01   1  96  .924      0.01  1  .924    
##  Carbon_t <= Carbon_p         0.07   1  96  .788      0.07  1  .788    
##  Carbon_t <= Tech             2.86   1  96  .094 .    2.86  1  .091 .  
##  Carbon_t <= ALL              2.70  10  96  .006 **  27.01 10  .003 ** 
##  ------------------------                                              
##  Tech <= iShares_GCE          0.69   1  96  .409      0.69  1  .407    
##  Tech <= Geop                 8.41   1  96  .005 **   8.41  1  .004 ** 
##  Tech <= EPU                  0.25   1  96  .617      0.25  1  .616    
##  Tech <= Energy_i             1.66   1  96  .201      1.66  1  .198    
##  Tech <= Energy_s             0.23   1  96  .631      0.23  1  .630    
##  Tech <= Green_e              2.22   1  96  .140      2.22  1  .137    
##  Tech <= Global_w             1.31   1  96  .255      1.31  1  .252    
##  Tech <= Natural_d            0.58   1  96  .449      0.58  1  .447    
##  Tech <= Carbon_p             2.09   1  96  .151      2.09  1  .148    
##  Tech <= Carbon_t             0.14   1  96  .714      0.14  1  .713    
##  Tech <= ALL                  2.72  10  96  .005 **  27.22 10  .002 ** 
## ───────────────────────────────────────────────────────────────────────
cat("\033[34mSP_GCE:", "\033[39m\n")
## [34mSP_GCE: [39m
vm = VAR(data_pred[1:110, c(16, 3:4, 6:13)], p = 1, type = 'both')
granger_causality(vm)
## 
## Granger Causality Test (Multivariate)
## 
## F test and Wald χ² test based on VAR(1) model:
## ────────────────────────────────────────────────────────────────────
##                               F df1 df2     p     Chisq df     p    
## ────────────────────────────────────────────────────────────────────
##  ---------------------                                              
##  SP_GCE <= Geop            0.03   1  96  .874      0.03  1  .874    
##  SP_GCE <= EPU             7.19   1  96  .009 **   7.19  1  .007 ** 
##  SP_GCE <= Energy_i       13.00   1  96 <.001 *** 13.00  1 <.001 ***
##  SP_GCE <= Energy_s        2.72   1  96  .102      2.72  1  .099 .  
##  SP_GCE <= Green_e         0.27   1  96  .603      0.27  1  .602    
##  SP_GCE <= Global_w        0.68   1  96  .411      0.68  1  .409    
##  SP_GCE <= Natural_d       0.34   1  96  .564      0.34  1  .563    
##  SP_GCE <= Carbon_p        0.51   1  96  .476      0.51  1  .474    
##  SP_GCE <= Carbon_t        0.18   1  96  .670      0.18  1  .669    
##  SP_GCE <= Tech            1.41   1  96  .238      1.41  1  .235    
##  SP_GCE <= ALL             3.29  10  96  .001 **  32.94 10 <.001 ***
##  ---------------------                                              
##  Geop <= SP_GCE            1.93   1  96  .168      1.93  1  .164    
##  Geop <= EPU               1.73   1  96  .192      1.73  1  .189    
##  Geop <= Energy_i          1.41   1  96  .238      1.41  1  .235    
##  Geop <= Energy_s          1.34   1  96  .250      1.34  1  .248    
##  Geop <= Green_e           0.19   1  96  .664      0.19  1  .663    
##  Geop <= Global_w          0.05   1  96  .825      0.05  1  .824    
##  Geop <= Natural_d         1.26   1  96  .264      1.26  1  .261    
##  Geop <= Carbon_p          1.08   1  96  .302      1.08  1  .300    
##  Geop <= Carbon_t          0.10   1  96  .747      0.10  1  .746    
##  Geop <= Tech              2.06   1  96  .154      2.06  1  .151    
##  Geop <= ALL               1.59  10  96  .122     15.88 10  .103    
##  ---------------------                                              
##  EPU <= SP_GCE             1.19   1  96  .279      1.19  1  .276    
##  EPU <= Geop               0.61   1  96  .435      0.61  1  .433    
##  EPU <= Energy_i           0.45   1  96  .503      0.45  1  .502    
##  EPU <= Energy_s           1.20   1  96  .275      1.20  1  .273    
##  EPU <= Green_e            1.96   1  96  .164      1.96  1  .161    
##  EPU <= Global_w           1.36   1  96  .247      1.36  1  .244    
##  EPU <= Natural_d          3.78   1  96  .055 .    3.78  1  .052 .  
##  EPU <= Carbon_p           1.24   1  96  .269      1.24  1  .266    
##  EPU <= Carbon_t           0.37   1  96  .542      0.37  1  .540    
##  EPU <= Tech               3.46   1  96  .066 .    3.46  1  .063 .  
##  EPU <= ALL                2.13  10  96  .029 *   21.27 10  .019 *  
##  ---------------------                                              
##  Energy_i <= SP_GCE        3.59   1  96  .061 .    3.59  1  .058 .  
##  Energy_i <= Geop          0.30   1  96  .583      0.30  1  .582    
##  Energy_i <= EPU           0.75   1  96  .390      0.75  1  .388    
##  Energy_i <= Energy_s     11.59   1  96 <.001 *** 11.59  1 <.001 ***
##  Energy_i <= Green_e       0.72   1  96  .398      0.72  1  .396    
##  Energy_i <= Global_w      0.16   1  96  .692      0.16  1  .691    
##  Energy_i <= Natural_d     0.13   1  96  .717      0.13  1  .716    
##  Energy_i <= Carbon_p      0.57   1  96  .452      0.57  1  .450    
##  Energy_i <= Carbon_t      1.58   1  96  .212      1.58  1  .209    
##  Energy_i <= Tech          2.77   1  96  .100 .    2.77  1  .096 .  
##  Energy_i <= ALL           2.24  10  96  .021 *   22.44 10  .013 *  
##  ---------------------                                              
##  Energy_s <= SP_GCE        0.55   1  96  .461      0.55  1  .459    
##  Energy_s <= Geop          0.51   1  96  .477      0.51  1  .476    
##  Energy_s <= EPU           0.03   1  96  .874      0.03  1  .874    
##  Energy_s <= Energy_i      0.43   1  96  .513      0.43  1  .512    
##  Energy_s <= Green_e       0.39   1  96  .532      0.39  1  .530    
##  Energy_s <= Global_w      0.66   1  96  .420      0.66  1  .418    
##  Energy_s <= Natural_d     0.35   1  96  .557      0.35  1  .556    
##  Energy_s <= Carbon_p      3.60   1  96  .061 .    3.60  1  .058 .  
##  Energy_s <= Carbon_t      0.49   1  96  .487      0.49  1  .485    
##  Energy_s <= Tech          3.53   1  96  .063 .    3.53  1  .060 .  
##  Energy_s <= ALL           1.54  10  96  .138     15.37 10  .119    
##  ---------------------                                              
##  Green_e <= SP_GCE         1.16   1  96  .285      1.16  1  .282    
##  Green_e <= Geop           2.05   1  96  .155      2.05  1  .152    
##  Green_e <= EPU            0.34   1  96  .561      0.34  1  .560    
##  Green_e <= Energy_i       0.00   1  96  .971      0.00  1  .971    
##  Green_e <= Energy_s       1.10   1  96  .297      1.10  1  .294    
##  Green_e <= Global_w       6.13   1  96  .015 *    6.13  1  .013 *  
##  Green_e <= Natural_d      1.01   1  96  .317      1.01  1  .314    
##  Green_e <= Carbon_p       1.92   1  96  .169      1.92  1  .166    
##  Green_e <= Carbon_t       0.35   1  96  .558      0.35  1  .556    
##  Green_e <= Tech           2.17   1  96  .144      2.17  1  .140    
##  Green_e <= ALL            1.90  10  96  .055 .   18.96 10  .041 *  
##  ---------------------                                              
##  Global_w <= SP_GCE        3.96   1  96  .049 *    3.96  1  .047 *  
##  Global_w <= Geop          1.74   1  96  .191      1.74  1  .187    
##  Global_w <= EPU           0.17   1  96  .678      0.17  1  .677    
##  Global_w <= Energy_i      5.65   1  96  .019 *    5.65  1  .017 *  
##  Global_w <= Energy_s      1.77   1  96  .186      1.77  1  .183    
##  Global_w <= Green_e       1.55   1  96  .217      1.55  1  .213    
##  Global_w <= Natural_d     0.00   1  96  .964      0.00  1  .964    
##  Global_w <= Carbon_p      0.00   1  96  .950      0.00  1  .950    
##  Global_w <= Carbon_t      2.90   1  96  .092 .    2.90  1  .089 .  
##  Global_w <= Tech          7.78   1  96  .006 **   7.78  1  .005 ** 
##  Global_w <= ALL           1.91  10  96  .052 .   19.13 10  .039 *  
##  ---------------------                                              
##  Natural_d <= SP_GCE       0.52   1  96  .471      0.52  1  .470    
##  Natural_d <= Geop         1.60   1  96  .209      1.60  1  .206    
##  Natural_d <= EPU          0.05   1  96  .824      0.05  1  .823    
##  Natural_d <= Energy_i     4.19   1  96  .043 *    4.19  1  .041 *  
##  Natural_d <= Energy_s     2.18   1  96  .144      2.18  1  .140    
##  Natural_d <= Green_e      2.67   1  96  .105      2.67  1  .102    
##  Natural_d <= Global_w     0.00   1  96  .969      0.00  1  .969    
##  Natural_d <= Carbon_p     0.11   1  96  .742      0.11  1  .741    
##  Natural_d <= Carbon_t     0.67   1  96  .416      0.67  1  .414    
##  Natural_d <= Tech         5.85   1  96  .017 *    5.85  1  .016 *  
##  Natural_d <= ALL          1.68  10  96  .096 .   16.81 10  .079 .  
##  ---------------------                                              
##  Carbon_p <= SP_GCE        1.94   1  96  .167      1.94  1  .164    
##  Carbon_p <= Geop          0.96   1  96  .330      0.96  1  .328    
##  Carbon_p <= EPU           1.60   1  96  .209      1.60  1  .206    
##  Carbon_p <= Energy_i      0.27   1  96  .604      0.27  1  .602    
##  Carbon_p <= Energy_s      0.16   1  96  .692      0.16  1  .692    
##  Carbon_p <= Green_e       0.54   1  96  .464      0.54  1  .462    
##  Carbon_p <= Global_w      0.36   1  96  .548      0.36  1  .546    
##  Carbon_p <= Natural_d     0.58   1  96  .449      0.58  1  .447    
##  Carbon_p <= Carbon_t      2.10   1  96  .150      2.10  1  .147    
##  Carbon_p <= Tech          4.68   1  96  .033 *    4.68  1  .031 *  
##  Carbon_p <= ALL           1.51  10  96  .147     15.11 10  .128    
##  ---------------------                                              
##  Carbon_t <= SP_GCE        0.24   1  96  .626      0.24  1  .625    
##  Carbon_t <= Geop          0.36   1  96  .552      0.36  1  .551    
##  Carbon_t <= EPU           0.46   1  96  .499      0.46  1  .498    
##  Carbon_t <= Energy_i      4.67   1  96  .033 *    4.67  1  .031 *  
##  Carbon_t <= Energy_s      3.20   1  96  .077 .    3.20  1  .074 .  
##  Carbon_t <= Green_e       0.49   1  96  .486      0.49  1  .484    
##  Carbon_t <= Global_w      6.16   1  96  .015 *    6.16  1  .013 *  
##  Carbon_t <= Natural_d     0.01   1  96  .917      0.01  1  .917    
##  Carbon_t <= Carbon_p      0.07   1  96  .797      0.07  1  .797    
##  Carbon_t <= Tech          2.82   1  96  .097 .    2.82  1  .093 .  
##  Carbon_t <= ALL           2.70  10  96  .006 **  26.97 10  .003 ** 
##  ---------------------                                              
##  Tech <= SP_GCE            0.67   1  96  .414      0.67  1  .412    
##  Tech <= Geop              8.39   1  96  .005 **   8.39  1  .004 ** 
##  Tech <= EPU               0.26   1  96  .608      0.26  1  .607    
##  Tech <= Energy_i          1.65   1  96  .202      1.65  1  .199    
##  Tech <= Energy_s          0.23   1  96  .631      0.23  1  .629    
##  Tech <= Green_e           2.21   1  96  .140      2.21  1  .137    
##  Tech <= Global_w          1.32   1  96  .253      1.32  1  .250    
##  Tech <= Natural_d         0.59   1  96  .445      0.59  1  .443    
##  Tech <= Carbon_p          2.14   1  96  .147      2.14  1  .143    
##  Tech <= Carbon_t          0.14   1  96  .714      0.14  1  .713    
##  Tech <= ALL               2.72  10  96  .006 **  27.20 10  .002 ** 
## ────────────────────────────────────────────────────────────────────
#ARIMA RENIXX
all_series <- ts(data_pred, start = c(2015, 1), frequency = 12)
for (i in 73:(sample - 1)) {
  arima <- auto.arima(all_series[1:i, c('Renixx')], seasonal = FALSE)
  pred_arima <- forecast::forecast(arima, h = 1)
  predictions_renixx[(i - 73 + 1), 1] <- pred_arima$mean
  lower_uni_renixx[(i - 73 + 1), 1] <- pred_arima$lower[, 2]
  up_uni_renixx[(i - 73 + 1), 1] <- pred_arima$upper[, 2]
}
#
#ARIMA MSCI_GAE
for (i in 73:(sample - 1)) {
  arima <- auto.arima(all_series[1:i, c('MSCI_GAE')], seasonal = FALSE)
  pred_arima <- forecast::forecast(arima, h = 1)
  predictions_MSCI_GAE[(i - 73 + 1), 1] <- pred_arima$mean
}
#
#ARIMA iShares_GCE
for (i in 73:(sample - 1)) {
  arima <- auto.arima(all_series[1:i, c('iShares_GCE')], seasonal = FALSE)
  pred_arima <- forecast::forecast(arima, h = 1)
  predictions_iShares_GCE[(i - 73 + 1), 1] <- pred_arima$mean
}
#
#ARIMA SP_GCE
for (i in 73:(sample - 1)) {
  arima <- auto.arima(all_series[1:i, c('SP_GCE')], seasonal = FALSE)
  pred_arima <- forecast::forecast(arima, h = 1)
  predictions_SP_GCE[(i - 73 + 1), 1] <- pred_arima$mean
}

This next chunk may take longer to execute than the previous ones, due to the computationally intensive estimation procedure.

#BSTS RENIXX
for (i in 73:(sample - 1)) {
  #Lag
  (ss <- list())
  ss <- AddAutoAr(ss, data_pred[1:i, c('Renixx')], lags = 1) #p=1 default
  ss <- AddLocalLevel(ss, data_pred[1:i, c('Renixx')])
  ss <- AddSeasonal(ss, data_pred[1:i, c('Renixx')], nseasons = 12)
  bsts.model <- bsts(
    data_pred[1:i, c('Renixx')],
    state.specification = ss,
    niter = 1000,
    seed = 009,
    #Number of MCMC draws
    ping = 0,
    expected.model.size = 3
  )
  burn <- SuggestBurn(0.1, bsts.model) #Suggested burning 0.1
  pred_bayesian <-
    predict(
      bsts.model,
      data_pred[(i + 1), ],
      horizon = 1,
      burn = burn,
      seed = 009,
      quantiles = c(.025, .975)
    )
  #Saving
  predictions_renixx[(i - 73 + 1), 2] <- pred_bayesian$mean
  lower_uni_renixx[(i - 73 + 1), 2] <- pred_bayesian$interval[1, ]
  up_uni_renixx[(i - 73 + 1), 2] <- pred_bayesian$interval[2, ]
}
#
#BSTS MSCI_GAE
for (i in 73:(sample - 1)) {
  #Lag
  (ss <- list())
  ss <- AddAutoAr(ss, data_pred[1:i, c('MSCI_GAE')], lags = 1) #p=1 default
  ss <- AddLocalLevel(ss, data_pred[1:i, c('MSCI_GAE')])
  ss <- AddSeasonal(ss, data_pred[1:i, c('MSCI_GAE')], nseasons = 12)
  bsts.model <- bsts(
    data_pred[1:i, c('MSCI_GAE')],
    state.specification = ss,
    niter = 1000,
    seed = 009,
    #Number of MCMC draws
    ping = 0,
    expected.model.size = 3
  )
  burn <- SuggestBurn(0.1, bsts.model) #Suggested burning 0.1
  pred_bayesian <-
    predict(
      bsts.model,
      data_pred[(i + 1), ],
      horizon = 1,
      burn = burn,
      seed = 009,
      quantiles = c(.025, .975)
    )
  #Saving
  predictions_MSCI_GAE[(i - 73 + 1), 2] <- pred_bayesian$mean
}
#
#BSTS iShares_GCE
for (i in 73:(sample - 1)) {
  #Lag
  (ss <- list())
  ss <- AddAutoAr(ss, data_pred[1:i, c('iShares_GCE')], lags = 1) #p=1 default
  ss <- AddLocalLevel(ss, data_pred[1:i, c('iShares_GCE')])
  ss <- AddSeasonal(ss, data_pred[1:i, c('iShares_GCE')], nseasons = 12)
  bsts.model <- bsts(
    data_pred[1:i, c('iShares_GCE')],
    state.specification = ss,
    niter = 1000,
    seed = 009,
    #Number of MCMC draws
    ping = 0,
    expected.model.size = 3
  )
  burn <- SuggestBurn(0.1, bsts.model) #Suggested burning 0.1
  pred_bayesian <-
    predict(
      bsts.model,
      data_pred[(i + 1), ],
      horizon = 1,
      burn = burn,
      seed = 009,
      quantiles = c(.025, .975)
    )
  #Saving
  predictions_iShares_GCE[(i - 73 + 1), 2] <- pred_bayesian$mean
}
#
#BSTS SP_GCE
for (i in 73:(sample - 1)) {
  #Lag
  (ss <- list())
  ss <- AddAutoAr(ss, data_pred[1:i, c('SP_GCE')], lags = 1) #p=1 default
  ss <- AddLocalLevel(ss, data_pred[1:i, c('SP_GCE')])
  ss <- AddSeasonal(ss, data_pred[1:i, c('SP_GCE')], nseasons = 12)
  bsts.model <- bsts(
    data_pred[1:i, c('SP_GCE')],
    state.specification = ss,
    niter = 1000,
    seed = 009,
    #Number of MCMC draws
    ping = 0,
    expected.model.size = 3
  )
  burn <- SuggestBurn(0.1, bsts.model) #Suggested burning 0.1
  pred_bayesian <-
    predict(
      bsts.model,
      data_pred[(i + 1), ],
      horizon = 1,
      burn = burn,
      seed = 009,
      quantiles = c(.025, .975)
    )
  #Saving
  predictions_SP_GCE[(i - 73 + 1), 2] <- pred_bayesian$mean
}

This next chunk may take longer to execute than the previous ones, due to the computationally intensive estimation procedure.

#ARIMAX with lagged values RENIXX
data_pred2 <-
  cbind(data_pred$Renixx, lag(data_pred[, c('Energy_i')]))[2:length(data_pred[, 1]), ] #Lagged values to prevent overconfidence
rownames(data_pred2) <- NULL
colnames(data_pred2) <- c('Renixx', 'Energy_i')
all_series <- ts(data_pred2, start = c(2015, 2), frequency = 12)
for (i in (73 - 1):(sample - 1 - 1)) {
  arimax <-
    auto.arima(all_series[1:i, 'Renixx'], seasonal = FALSE, xreg = all_series[1:i, c('Energy_i'), drop = FALSE])
  new <- matrix(all_series[(i + 1), c('Energy_i')], ncol = 1)
  colnames(new) <- c('Energy_i')
  pred_arimax <- forecast::forecast(arimax, xreg = new)
  predictions_renixx[(i - (73 - 1) + 1), 3] <- pred_arimax$mean
  lower_uni_renixx[(i - (73 - 1) + 1), 3] <- pred_arimax$lower[, 2]
  up_uni_renixx[(i - (73 - 1) + 1), 3] <- pred_arimax$upper[, 2]
}
#
#ARIMAX with lagged values MSCI_GAE
data_pred2 <-
  cbind(data_pred$MSCI_GAE, lag(data_pred[, 'Energy_i']))[2:length(data_pred[, 1]), ] #Lagged values to prevent overconfidence
rownames(data_pred2) <- NULL
colnames(data_pred2) <- c('MSCI_GAE', 'Energy_i')
all_series <- ts(data_pred2, start = c(2015, 2), frequency = 12)
for (i in (73 - 1):(sample - 1 - 1)) {
  arimax <-
    auto.arima(all_series[1:i, 'MSCI_GAE'], seasonal = FALSE, xreg = all_series[1:i, c('Energy_i'), drop = FALSE])
  new <- matrix(all_series[(i + 1), c('Energy_i')], ncol = 1)
  colnames(new) <- c('Energy_i')
  pred_arimax <- forecast::forecast(arimax, xreg = new)
  predictions_MSCI_GAE[(i - (73 - 1) + 1), 3] <- pred_arimax$mean
  lower_uni_MSCI_GAE[(i - (73 - 1) + 1), 3] <- pred_arimax$lower[, 2]
  up_uni_MSCI_GAE[(i - (73 - 1) + 1), 3] <- pred_arimax$upper[, 2]
}
#
#ARIMAX with lagged values iShares_GCE
data_pred2 <-
  cbind(data_pred$iShares_GCE, lag(data_pred[, 'Energy_i']))[2:length(data_pred[, 1]), ] #Lagged values to prevent overconfidence
rownames(data_pred2) <- NULL
colnames(data_pred2) <- c('iShares_GCE', 'Energy_i')
all_series <- ts(data_pred2, start = c(2015, 2), frequency = 12)
for (i in (73 - 1):(sample - 1 - 1)) {
  arimax <-
    auto.arima(all_series[1:i, 'iShares_GCE'], seasonal = FALSE, xreg = all_series[1:i, c('Energy_i'), drop = FALSE])
  new <- matrix(all_series[(i + 1), c('Energy_i')], ncol = 1)
  colnames(new) <- c('Energy_i')
  pred_arimax <- forecast::forecast(arimax, xreg = new)
  predictions_iShares_GCE[(i - (73 - 1) + 1), 3] <- pred_arimax$mean
  lower_uni_iShares_GCE[(i - (73 - 1) + 1), 3] <- pred_arimax$lower[, 2]
  up_uni_iShares_GCE[(i - (73 - 1) + 1), 3] <- pred_arimax$upper[, 2]
}
#
#ARIMAX with lagged values SP_GCE
data_pred2 <-
  cbind(data_pred$SP_GCE, lag(data_pred[, 'Energy_i']))[2:length(data_pred[, 1]), ] #Lagged values to prevent overconfidence
rownames(data_pred2) <- NULL
colnames(data_pred2) <- c('SP_GCE', 'Energy_i')
all_series <- ts(data_pred2, start = c(2015, 2), frequency = 12)
for (i in (73 - 1):(sample - 1 - 1)) {
  arimax <-
    auto.arima(all_series[1:i, 'SP_GCE'], seasonal = FALSE, xreg = all_series[1:i, c('Energy_i'), drop = FALSE])
  new <- matrix(all_series[(i + 1), c('Energy_i')], ncol = 1)
  colnames(new) <- c('Energy_i')
  pred_arimax <- forecast::forecast(arimax, xreg = new)
  predictions_SP_GCE[(i - (73 - 1) + 1), 3] <- pred_arimax$mean
  lower_uni_SP_GCE[(i - (73 - 1) + 1), 3] <- pred_arimax$lower[, 2]
  up_uni_SP_GCE[(i - (73 - 1) + 1), 3] <- pred_arimax$upper[, 2]
}

This next chunk may take longer to execute than the previous ones, due to the computationally intensive estimation procedure.

#ARIMAX with current values RENIXX
all_series <- ts(data_pred, start = c(2015, 1), frequency = 12)
for (i in 73:(sample - 1)) {
  x <- matrix(all_series[1:i, c('Energy_i')], ncol =
                1)
  colnames(x) <- c('Energy_i')
  arimax <- auto.arima(all_series[1:i, c('Renixx')], xreg = x, seasonal = FALSE)
  new <-
    matrix(all_series[(i + 1), c('Energy_i')], ncol =
             1)
  colnames(new) <- c('Energy_i')
  pred_arimax <- forecast::forecast(arimax, xreg = new)
  predictions_renixx[(i - 73 + 1), 4] <- pred_arimax$mean
  lower_uni_renixx[(i - 73 + 1), 4] <- pred_arimax$lower[, 2]
  up_uni_renixx[(i - 73 + 1), 4] <- pred_arimax$upper[, 2]
}
#
#ARIMAX with current values MSCI_GAE
for (i in 73:(sample - 1)) {
  x <- matrix(all_series[1:i, c('Energy_i')], ncol =
                1)
  colnames(x) <- c('Energy_i')
  arimax <- auto.arima(all_series[1:i, c('MSCI_GAE')], xreg = x, seasonal = FALSE)
  new <-
    matrix(all_series[(i + 1), c('Energy_i')], ncol =
             1)
  colnames(new) <- c('Energy_i')
  pred_arimax <- forecast::forecast(arimax, xreg = new)
  predictions_MSCI_GAE[(i - 73 + 1), 4] <- pred_arimax$mean
  lower_uni_MSCI_GAE[(i - 73 + 1), 4] <- pred_arimax$lower[, 2]
  up_uni_MSCI_GAE[(i - 73 + 1), 4] <- pred_arimax$upper[, 2]
}
#
#ARIMAX with current values iShares_GCE
for (i in 73:(sample - 1)) {
  x <- matrix(all_series[1:i, c('Energy_i')], ncol =
                1)
  colnames(x) <- c('Energy_i')
  arimax <- auto.arima(all_series[1:i, c('iShares_GCE')], xreg = x, seasonal = FALSE)
  new <-
    matrix(all_series[(i + 1), c('Energy_i')], ncol =
             1)
  colnames(new) <- c('Energy_i')
  pred_arimax <- forecast::forecast(arimax, xreg = new)
  predictions_iShares_GCE[(i - 73 + 1), 4] <- pred_arimax$mean
  lower_uni_iShares_GCE[(i - 73 + 1), 4] <- pred_arimax$lower[, 2]
  up_uni_iShares_GCE[(i - 73 + 1), 4] <- pred_arimax$upper[, 2]
}
#
#ARIMAX with current values SP_GCE
for (i in 73:(sample - 1)) {
  x <- matrix(all_series[1:i, c('Energy_i')], ncol =
                1)
  colnames(x) <- c('Energy_i')
  arimax <- auto.arima(all_series[1:i, c('SP_GCE')], xreg = x, seasonal = FALSE)
  new <-
    matrix(all_series[(i + 1), c('Energy_i')], ncol =
             1)
  colnames(new) <- c('Energy_i')
  pred_arimax <- forecast::forecast(arimax, xreg = new)
  predictions_SP_GCE[(i - 73 + 1), 4] <- pred_arimax$mean
  lower_uni_SP_GCE[(i - 73 + 1), 4] <- pred_arimax$lower[, 2]
  up_uni_SP_GCE[(i - 73 + 1), 4] <- pred_arimax$upper[, 2]
}

Below, the metrics of the model are reported, along with some robustness checks (Table 6 in the article).

#Vector of predictions
predictions_uni <- rbind(
  predictions_renixx,
  predictions_MSCI_GAE,
  predictions_iShares_GCE,
  predictions_SP_GCE
)
#Vector of real data
real_data <- c(data_pred$Renixx[74:110],
               data_pred$MSCI_GAE[74:110],
               data_pred$iShares_GCE[74:110],
               data_pred$SP_GCE[74:110])
#SMAPE
cat("SMAPE ARIMA:",
    Metrics::smape(predictions_uni[, 1], real_data) * 100,
    "\n") #By 100 to obtain the percentage
## SMAPE ARIMA: 1.4989211
cat("SMAPE BSTS:",
    Metrics::smape(predictions_uni[, 2], real_data) * 100,
    "\n")
## SMAPE BSTS: 1.36563165
cat("SMAPE ARIMAX_t1:",
    Metrics::smape(predictions_uni[, 3], real_data) * 100,
    "\n")
## SMAPE ARIMAX_t1: 1.31293627
cat("SMAPE ARIMAX_t:",
    Metrics::smape(predictions_uni[, 4], real_data) * 100,
    "\n")
## SMAPE ARIMAX_t: 1.26852154
#RMSE
cat("RMSE ARIMA:",
    Metrics::rmse(predictions_uni[, 1], real_data),
    "\n")
## RMSE ARIMA: 0.0905562743
cat("RMSE BSTS:", Metrics::rmse(predictions_uni[, 2], real_data), "\n")
## RMSE BSTS: 0.0814020975
cat("RMSE ARIMAX_t1:",
    Metrics::rmse(predictions_uni[, 3], real_data),
    "\n")
## RMSE ARIMAX_t1: 0.0771483961
cat("RMSE ARIMAX_t:",
    Metrics::rmse(predictions_uni[, 4], real_data),
    "\n")
## RMSE ARIMAX_t: 0.0728225047
set.seed(005)
#DM against ARIMAX with actual values (t)
cat("\033[34mP-value of the Diebold-Mariano test: ARIMA",
    "\033[39m\n")
## [34mP-value of the Diebold-Mariano test: ARIMA [39m
DM.test(
  predictions_uni[, 4],
  predictions_uni[, 1],
  real_data,
  loss.type = "SE",
  1,
  c = FALSE,
  H1 = "more"
)
## 
##  Diebold-Mariano test
## 
## data:  predictions_uni[, 4] and predictions_uni[, 1] and real_data
## statistic = -2.774686, forecast horizon = 1, p-value = 0.00276275
## alternative hypothesis: Forecast f1 is more accurate than f2.
cat("\033[34mP-value of the Diebold-Mariano test: BSTS",
    "\033[39m\n")
## [34mP-value of the Diebold-Mariano test: BSTS [39m
DM.test(
  predictions_uni[, 4],
  predictions_uni[, 2],
  real_data,
  loss.type = "SE",
  1,
  c = FALSE,
  H1 = "more"
)
## 
##  Diebold-Mariano test
## 
## data:  predictions_uni[, 4] and predictions_uni[, 2] and real_data
## statistic = -2.299387, forecast horizon = 1, p-value = 0.0107415
## alternative hypothesis: Forecast f1 is more accurate than f2.
cat("\033[34mP-value of the Diebold-Mariano test: ARIMAX_t1",
    "\033[39m\n")
## [34mP-value of the Diebold-Mariano test: ARIMAX_t1 [39m
DM.test(
  predictions_uni[, 4],
  predictions_uni[, 3],
  real_data,
  loss.type = "SE",
  1,
  c = FALSE,
  H1 = "more"
)
## 
##  Diebold-Mariano test
## 
## data:  predictions_uni[, 4] and predictions_uni[, 3] and real_data
## statistic = -1.009512, forecast horizon = 1, p-value = 0.156365
## alternative hypothesis: Forecast f1 is more accurate than f2.
#MCS
cat("\033[34mModel Confidence Set procedure", "\033[39m\n")
## [34mModel Confidence Set procedure [39m
Loss <- as.matrix((predictions_uni[, c(1, 2, 3, 4)] - real_data)^2)
suppressWarnings({
  MCSprocedure(Loss,
               alpha = 0.15,
               B = 5000,
               statistic = "TR")
})
## 
## Model BSTS eliminated 2025-12-10 12:29:18.53121
## Model ARIMA eliminated 2025-12-10 12:29:20.833602
## ###########################################################################################################################
## Superior Set Model created   :
##           Rank_M         v_M  MCS_M Rank_R         v_R  MCS_R          Loss
## ARIMAX_t1      2  1.11113149 0.3196      2  1.11113149 0.3196 0.00595187503
## ARIMAX_t       1 -1.11113149 1.0000      1 -1.11113149 1.0000 0.00530311719
## p-value  :
## [1] 0.3196
## 
## ###########################################################################################################################
## 
## ------------------------------------------
## -          Superior Set of Models        -
## ------------------------------------------
##           Rank_M         v_M  MCS_M Rank_R         v_R  MCS_R          Loss
## ARIMAX_t1      2  1.11113149 0.3196      2  1.11113149 0.3196 0.00595187503
## ARIMAX_t       1 -1.11113149 1.0000      1 -1.11113149 1.0000 0.00530311719
## 
## Details
## ------------------------------------------
## 
## Number of eliminated models  :   2
## Statistic    :   TR
## Elapsed Time :   Time difference of 6.28184056 secs

The ARIMAX models, with the energy index as an exogenous variable, are identified as the best models based on both RMSE and MAPE, as well as robustness checks such as the Diebold-Mariano test and the Model Confidence Set procedure. Finally, the results of the best model are presented below (Figure 5 in the article).

#Plot RENIXX
p1 <- ggplot(data_pred[65:sample, ], aes(x = as.Date(Date), y = Renixx), ) +
  geom_line(aes(
    y = c(data_pred$Renixx[65:73], predictions_renixx[, 3]),
    color = 'Prediction'
  )) +
  geom_ribbon(
    aes(
      ymin = c(data_pred$Renixx[65:73], lower_uni_renixx[, 3]),
      ymax = c(data_pred$Renixx[65:73], up_uni_renixx[, 3])
    ),
    alpha = 0.1,
    fill = "green",
    color = "green",
    linetype = "dotted"
  ) +
  geom_line(aes(color = 'Real data')) +
  labs(title = 'RENIXX Index') + xlab('Time') + ylab('Log-price') +
  scale_color_manual(
    name = 'Legend',
    breaks = c('Real data', 'Prediction'),
    values = c(
      'Real data' = 'black',
      'Prediction' = 'darkgreen'
    )
  )
#Plot MSCI_GAE
p2 <- ggplot(data_pred[65:sample, ], aes(x = as.Date(Date), y = MSCI_GAE), ) +
  geom_line(aes(
    y = c(data_pred$MSCI_GAE[65:73], predictions_MSCI_GAE[, 3]),
    color = 'Prediction'
  )) +
  geom_ribbon(
    aes(
      ymin = c(data_pred$MSCI_GAE[65:73], lower_uni_MSCI_GAE[, 3]),
      ymax = c(data_pred$MSCI_GAE[65:73], up_uni_MSCI_GAE[, 3])
    ),
    alpha = 0.1,
    fill = "blue",
    color = "blue",
    linetype = "dotted"
  ) +
  geom_line(aes(color = 'Real data')) +
  labs(title = 'MSCI Global Alternative Energy Index') + xlab('Time') + ylab('Log-price') +
  scale_color_manual(
    name = 'Legend',
    breaks = c('Real data', 'Prediction'),
    values = c(
      'Real data' = 'black',
      'Prediction' = 'darkblue'
    )
  )
#Plot iShares_GCE
p3 <- ggplot(data_pred[65:sample, ], aes(x = as.Date(Date), y = iShares_GCE), ) +
  geom_line(aes(
    y = c(data_pred$iShares_GCE[65:73], predictions_iShares_GCE[, 3]),
    color = 'Prediction'
  )) +
  geom_ribbon(
    aes(
      ymin = c(data_pred$iShares_GCE[65:73], lower_uni_iShares_GCE[, 3]),
      ymax = c(data_pred$iShares_GCE[65:73], up_uni_iShares_GCE[, 3])
    ),
    alpha = 0.1,
    fill = "red",
    color = "red",
    linetype = "dotted"
  ) +
  geom_line(aes(color = 'Real data')) +
  labs(title = 'IShares Global Clean Energy Index') + xlab('Time') + ylab('Log-price') +
  scale_color_manual(
    name = 'Legend',
    breaks = c('Real data', 'Prediction'),
    values = c(
      'Real data' = 'black',
      'Prediction' = 'darkred'
    )
  )
#Plot SP_GCE
p4 <- ggplot(data_pred[65:sample, ], aes(x = as.Date(Date), y = SP_GCE), ) +
  geom_line(aes(
    y = c(data_pred$SP_GCE[65:73], predictions_SP_GCE[, 3]),
    color = 'Prediction'
  )) +
  geom_ribbon(
    aes(
      ymin = c(data_pred$SP_GCE[65:73], lower_uni_SP_GCE[, 3]),
      ymax = c(data_pred$SP_GCE[65:73], up_uni_SP_GCE[, 3])
    ),
    alpha = 0.1,
    fill = "violet",
    color = "violet",
    linetype = "dotted"
  ) +
  geom_line(aes(color = 'Real data')) +
  labs(title = 'S&P Global Clean Energy Index') + xlab('Time') + ylab('Log-price') +
  scale_color_manual(
    name = 'Legend',
    breaks = c('Real data', 'Prediction'),
    values = c('Real data' = 'black', 'Prediction' = 'purple')
  )
#Overall plot
grid.arrange(p1, p2, p3, p4, ncol = 2, top = textGrob("", gp = gpar(fontsize = 20, font = 3)))

Below are the summary statistics for the considered dataset (Table 1 in the article).

#Dataset
data_pred <-
  data.frame(
    renixx_m$close,
    MSCI_m$Close,
    geo_m$GPRH,
    GEPUCURRENT$GEPUCURRENT,
    Oil_m$Close,
    Fs_m$OFR.FSI,
    Energy_i[, 3],
    Energy_s[, 3],
    green_e[, 3],
    global_w[, 3],
    Natural[, 3],
    Carbon_p[, 3],
    Tax[, 3],
    Tech_d[, 3],
    green_e[, 1] #Date
  )[1:230, ] #2005/01/01 2024/02/01
colnames(data_pred) <-
  c(
    'Renixx',
    'MSCI',
    'Geop',
    'EPU',
    'Oil_p',
    'FSI',
    'Energy_i',
    'Energy_s',
    'Green_e',
    'Global_w',
    'Natural_d',
    'Carbon_p',
    'Carbon_t',
    'Tech',
    'Date'
  )
rownames(data_pred) <- NULL
options(digits = 9) #For 3-digits results
summary(data_pred[, c(1, 5, 2, 4, 6, 3, 7:14)])
##      Renixx             Oil_p               MSCI               EPU          
##  Min.   : 154.172   Min.   : 24.5024   Min.   : 585.702   Min.   : 49.2283  
##  1st Qu.: 434.340   1st Qu.: 49.0873   1st Qu.: 982.644   1st Qu.:102.5598  
##  Median : 575.251   Median : 57.9886   Median :1291.398   Median :145.2489  
##  Mean   : 776.210   Mean   : 61.5484   Mean   :1496.991   Mean   :162.3685  
##  3rd Qu.:1102.412   3rd Qu.: 77.5155   3rd Qu.:1865.287   3rd Qu.:220.4218  
##  Max.   :2169.222   Max.   :111.2582   Max.   :3042.451   Max.   :431.7576  
##       FSI                 Geop             Energy_i          Energy_s       
##  Min.   :-4.994368   Min.   : 56.9939   Min.   : 7391.0   Min.   : 1645.00  
##  1st Qu.:-3.127738   1st Qu.: 81.0432   1st Qu.:11086.0   1st Qu.: 6170.00  
##  Median :-1.620596   Median : 90.4331   Median :13857.5   Median : 7404.00  
##  Mean   :-0.322391   Mean   : 96.5832   Mean   :16013.9   Mean   : 8558.97  
##  3rd Qu.: 0.771736   3rd Qu.:103.8940   3rd Qu.:20232.5   3rd Qu.: 9872.00  
##  Max.   :24.748391   Max.   :306.9907   Max.   :36955.0   Max.   :20567.00  
##     Green_e          Global_w         Natural_d         Carbon_p      
##  Min.   : 98425   Min.   : 148245   Min.   : 58999   Min.   :15210.0  
##  1st Qu.:156751   1st Qu.: 321197   1st Qu.: 93416   1st Qu.:27736.0  
##  Median :185914   Median : 395320   Median :113082   Median :35789.0  
##  Mean   :194267   Mean   : 579230   Mean   :121398   Mean   :40569.8  
##  3rd Qu.:229659   3rd Qu.: 691810   3rd Qu.:132749   3rd Qu.:55473.0  
##  Max.   :364538   Max.   :2470750   Max.   :491664   Max.   :89473.0  
##     Carbon_t             Tech       
##  Min.   :  6515.0   Min.   :156224  
##  1st Qu.: 26063.0   1st Qu.:230225  
##  Median : 32579.0   Median :271336  
##  Mean   : 43981.8   Mean   :296504  
##  3rd Qu.: 52127.0   3rd Qu.:335059  
##  Max.   :325795.0   Max.   :666008

To ensure full reproducibility, the following chunk displays the R version, platform details, and package versions used in this analysis.

sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 20.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3;  LAPACK version 3.9.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## time zone: UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] parallel  grid      stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] vars_1.6-1          lmtest_0.9-40       strucchange_1.5-4  
##  [4] sandwich_3.1-1      urca_1.3-4          tseries_0.10-58    
##  [7] rugarch_1.5-4       readxl_1.4.5        multDM_1.1.5       
## [10] Metrics_0.1.4       MCS_0.1.3           lubridate_1.9.4    
## [13] gridExtra_2.3       ghyp_1.6.5          MASS_7.3-65        
## [16] numDeriv_2016.8-1.1 ggfortify_0.4.19    forecast_8.24.0    
## [19] fDMA_2.2.8          exuber_1.1.0        cpm_2.3            
## [22] bsts_0.9.11         xts_0.14.1          zoo_1.8-14         
## [25] BoomSpikeSlab_1.2.7 Boom_0.9.16         cowplot_1.2.0      
## [28] ggplot2_4.0.0       interactions_1.2.0  lmerTest_3.1-3     
## [31] lme4_1.1-38         Matrix_1.7-4        performance_0.15.3 
## [34] effectsize_1.0.1    emmeans_2.0.0       data.table_1.17.8  
## [37] forcats_1.0.1       stringr_1.6.0       tidyr_1.3.1        
## [40] dplyr_1.1.4         bruceR_2025.8       bbdetection_1.0    
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3          rstudioapi_0.17.1          
##   [3] jsonlite_2.0.0              cvar_0.5                   
##   [5] datawizard_1.3.0            magrittr_2.0.4             
##   [7] estimability_1.5.1          farver_2.1.2               
##   [9] nloptr_2.2.1                rmarkdown_2.30             
##  [11] vctrs_0.6.5                 minqa_1.2.8                
##  [13] htmltools_0.5.9             progress_1.2.3             
##  [15] itertools_0.1-3             curl_7.0.0                 
##  [17] truncnorm_1.0-9             broom_1.0.11               
##  [19] cellranger_1.1.0            TTR_0.24.4                 
##  [21] pracma_2.4.6                sass_0.4.10                
##  [23] parallelly_1.45.1           KernSmooth_2.23-26         
##  [25] bslib_0.9.0                 cachem_1.1.0               
##  [27] lifecycle_1.0.4             iterators_1.0.14           
##  [29] pkgconfig_2.0.3             R6_2.6.1                   
##  [31] fastmap_1.2.0               rbibutils_2.4              
##  [33] future_1.68.0               digest_0.6.39              
##  [35] colorspace_2.1-2            spatial_7.3-18             
##  [37] furrr_0.3.1                 GeneralizedHyperbolic_0.8-7
##  [39] labeling_0.4.3              timechange_0.3.0           
##  [41] compiler_4.5.2              rngtools_1.5.2             
##  [43] withr_3.0.2                 doParallel_1.0.17          
##  [45] pander_0.6.6                S7_0.2.1                   
##  [47] backports_1.5.0             psych_2.5.6                
##  [49] gplots_3.3.0                broom.mixed_0.2.9.6        
##  [51] fBasics_4041.97             gtools_3.9.5               
##  [53] caTools_1.18.3              tools_4.5.2                
##  [55] quantmod_0.4.28             future.apply_1.20.0        
##  [57] nnet_7.3-20                 glue_1.8.0                 
##  [59] quadprog_1.5-8              nlme_3.1-168               
##  [61] fGarch_4033.92              generics_0.1.4             
##  [63] snow_0.4-4                  gtable_0.3.6               
##  [65] hms_1.1.4                   foreach_1.5.2              
##  [67] pillar_1.11.1               splines_4.5.2              
##  [69] lattice_0.22-7              ks_1.15.1                  
##  [71] tidyselect_1.2.1            SkewHyperbolic_0.4-2       
##  [73] knitr_1.50                  reformulas_0.4.2           
##  [75] xfun_0.54                   timeDate_4051.111          
##  [77] Rsolnp_2.0.1                MTS_1.2.1                  
##  [79] stringi_1.8.7               yaml_2.3.11                
##  [81] boot_1.3-32                 DistributionUtils_0.6-2    
##  [83] evaluate_1.0.5              codetools_0.2-20           
##  [85] timeSeries_4041.111         tibble_3.3.0               
##  [87] cli_3.6.5                   spd_2.0-1                  
##  [89] xtable_1.8-4                parameters_0.28.3          
##  [91] Rdpack_2.6.4                jquerylib_0.1.4            
##  [93] Rcpp_1.1.0                  gbutils_0.5                
##  [95] doSNOW_1.0.20               globals_0.18.0             
##  [97] coda_0.19-4.1               png_0.1-8                  
##  [99] fracdiff_1.5-3              prettyunits_1.2.0          
## [101] mclust_6.1.2                bayestestR_0.17.0          
## [103] doRNG_1.8.6.2               bitops_1.0-9               
## [105] listenv_0.10.0              mvtnorm_1.3-3              
## [107] scales_1.4.0                insight_1.4.4              
## [109] purrr_1.2.0                 crayon_1.5.3               
## [111] rlang_1.1.6                 mnormt_2.1.1               
## [113] jtools_2.3.0